HCP-D Tract Profiles
Audrey Luo
cohortfile <- read.csv("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/cohortfiles/dti_fa/dti_fa_cohortfile.csv")
all_subjects <- fread("/cbica/projects/luo_wm_dev/input/HCPD/HCPD_tractprofiles/all_subjects/collated_tract_profiles.tsv")
all_subjects <- all_subjects %>%
mutate(
hemi = ifelse(grepl("_L", tractID), "Left", "Right"),
tractID = gsub("_[LR]", "", tractID),
tractID = case_when(
tractID == "ATR" ~ "Anterior Thalamic Radiation",
tractID == "CGC" ~ "Cingulum Cingulate",
tractID == "CST" ~ "Corticospinal Tract",
tractID == "IFO" ~ "Inferior Fronto-occipital Fasciculus",
tractID == "ILF" ~ "Inferior Longitudinal Fasciculus",
tractID == "SLF" ~ "Superior Longitudinal Fasciculus",
tractID == "ARC" ~ "Arcuate Fasciculus",
tractID == "UNC" ~ "Uncinate Fasciculus",
tractID == "FA" ~ "Forceps Minor",
tractID == "FP" ~ "Forceps Major",
tractID == "pARC" ~ "Posterior Arcuate",
tractID == "VOF" ~ "Vertical Occipital Fasciculus",
TRUE ~ tractID
),
tract_hemi = paste(hemi, tractID),
tract_hemi = gsub("Right Forceps", "Forceps", tract_hemi),
tract_hemi = str_replace_all(tract_hemi, " ", "_")
)
all_subjects$hemi[all_subjects$tract_hemi == "Forceps_Minor" | all_subjects$tract_hemi == "Forceps_Major"] <- NA
all_subjects$subjectID <- as.factor(all_subjects$subjectID)
cohortfile <- cohortfile %>% rename(subjectID=sub) %>% select(subjectID, age, sex, mean_fd)Tract Names and Orientations
Note that I used an older version of pyAFQ here, so I had to flip some tracts to ensure consistent orientations. I will rerun tract profiles using most updated pyAFQ once newest qsiprep is released.
| Abbrev | Full Name | Node 0 | Node 99 | Need to flip? | Final Orientation |
|---|---|---|---|---|---|
| ATR | anterior thalamic radiation | A (Frontal) | P (Thalamus) | 0 | A - P (Frontal - Thalamus) |
| CGC | cingulum cinguluate | P | A | 1 | A - P |
| IFO | inferior fronto-occipital fasciculus | P (Occipital) | A (Frontal) | 1 | A - P (Frontal - Occipital) |
| ILF | inferior longitudinal fasciculus | P (Occipital) | A (Temporal) | 1 | A - P (Temporal - Occipital) |
| SLF | superior longitudinal fasciculus | A | P | 0 | A - P |
| ARC | arcuate | A (Frontal) | P (Temporal) | 0 | A - P (Frontal - Temporal) |
| UNC | uncinate fasciculus | P (Temporal) | A (Frontal) | 1 | A - P (Frontal - Temporal) |
| FA | forceps minor | R | L | 0 | R - L |
| FP | forceps major | R | L | 0 | R - L |
| pARC | posterior arcuate | S | I | 0 | S - I |
| CST | corticospinal | I (Brainstem) | S (Motor cortex) | 1 | S - I (Motor - Brainstem) |
| VOF | vertical occipital fasciculus | S | I | 0 | S - I |
# fix tract orientations
tract_profiles <- all_subjects %>%
group_by(tractID) %>%
mutate(nodeID = ifelse(tractID %in% c("Cingulum Cingulate", "Inferior Fronto-occipital Fasciculus", "Inferior Longitudinal Fasciculus", "Uncinate Fasciculus", "Corticospinal Tract"), max(nodeID) - nodeID, nodeID))
# label main orientation
tract_profiles <- tract_profiles %>%
mutate(main_orientation = case_when(
tractID %in% c("Anterior Thalamic Radiation", "Cingulum Cingulate", "Inferior Fronto-occipital Fasciculus",
"Inferior Longitudinal Fasciculus", "Superior Longitudinal Fasciculus") ~ "AP",
tractID %in% c("Arcuate Fasciculus", "Uncinate Fasciculus") ~ "AP_frontal_temporal",
tractID %in% c("Forceps Minor", "Forceps Major") ~ "RL",
tractID %in% c("Corticospinal Tract", "Posterior Arcuate","Vertical Occipital Fasciculus") ~ "SI",
TRUE ~ NA_character_
))
tract_profiles$main_orientation <- as.factor(tract_profiles$main_orientation)
#write.table(tract_profiles, "/cbica/projects/luo_wm_dev/input/HCPD/HCPD_tractprofiles/all_subjects/collated_tract_profiles_reoriented.tsv")1. Tract profiles
# Func to visualize GAM model outputs
# adapted from https://github.com/bart-larsen/GAMM-Tutorial/blob/master/GAMM_tutorial.Rmd
visualize_model_hemis <- function(modobj, modobj_rh = NULL, both_hemis = F, smooth_var, int_var, group_var, plabels = NULL,check_diagnostics = F,derivative_plot = F, custom_color, custom_color_rh = NULL) { # modobj corresponds to lh if plotting both hemispheres
this_font_size = font_size*1.25
if (both_hemis == T) { # if want to plot model for left and right tract
if (is.null(modobj_rh)) {
stop("Please provide models for both hemispheres")
}
if (any(class(modobj)=="gam")) {
model_lh <- modobj
model_rh <- modobj_rh
} else if (class(modobj$gam)=="gam") {
model_lh <- modobj$gam
model_rh <- modobj_rh$gam
} else {
stop("Can't find a gam object to plot")
}
s_lh <- summary(model_lh)
s_rh <- summary(model_rh)
## Generate custom line plot
np <- 100 #number of predicted values - originally 100000 but we're doing only 100 discrete nodes...
df_lh <- model_lh$model
df_lh <- df_lh %>% mutate(hemi = "Left")
df_rh <- model_rh$model
df_rh <- df_rh %>% mutate(hemi = "Right")
theseVars <- attr(model_lh$terms,"term.labels") # vars are identical for lh and rh
varClasses <- attr(model_lh$terms,"dataClasses")
thisResp <- as.character(model_lh$terms[[2]])
# No interaction variable, just produce a single line plot
thisPred_lh <- data.frame(init = rep(0,np))
thisPred_rh <- data.frame(init = rep(0,np))
# lh
for (v in c(1:length(theseVars))) {
thisVar_lh <- theseVars[[v]]
thisClass_lh <- varClasses[thisVar_lh]
if (thisVar_lh == smooth_var) {
thisPred_lh[,smooth_var] = seq(min(df_lh[,smooth_var],na.rm = T),max(df_lh[,smooth_var],na.rm = T), length.out = np)
} else {
switch (thisClass_lh,
"numeric" = {thisPred_lh[,thisVar] = median(df_lh[,thisVar])},
"factor" = {thisPred_lh[,thisVar] = levels(df_lh[,thisVar])[[1]]},
"ordered" = {thisPred_lh[,thisVar] = levels(df_lh[,thisVar])[[1]]}
)
}
}
pred_lh <- thisPred_lh %>% select(-init)
p_lh <-data.frame(predict(model_lh,pred_lh,se.fit = T))
pred_lh <- cbind(pred_lh,p_lh)
pred_lh$selo <- pred_lh$fit - 2*pred_lh$se.fit
pred_lh$sehi <- pred_lh$fit + 2*pred_lh$se.fit
pred_lh[,group_var] = NA
pred_lh[,thisResp] = 1
pred_lh <- pred_lh %>% mutate(hemi = "Left")
# rh
for (v in c(1:length(theseVars))) {
thisVar_rh <- theseVars[[v]]
thisClass_rh <- varClasses[thisVar_rh]
if (thisVar_rh == smooth_var) {
thisPred_rh[,smooth_var] = seq(min(df_rh[,smooth_var],na.rm = T),max(df_rh[,smooth_var],na.rm = T), length.out = np)
} else {
switch (thisClass_rh,
"numeric" = {thisPred_rh[,thisVar] = median(df_rh[,thisVar])},
"factor" = {thisPred_rh[,thisVar] = levels(df_rh[,thisVar])[[1]]},
"ordered" = {thisPred_rh[,thisVar] = levels(df_rh)[[1]]}
)
}
}
pred_rh <- thisPred_rh %>% select(-init)
p_rh <-data.frame(predict(model_rh,pred_rh,se.fit = T))
pred_rh <- cbind(pred_rh,p_rh)
pred_rh$selo <- pred_rh$fit - 2*pred_rh$se.fit
pred_rh$sehi <- pred_rh$fit + 2*pred_rh$se.fit
pred_rh[,group_var] = NA
pred_rh[,thisResp] = 1
pred_rh <- pred_rh %>% mutate(hemi = "Right")
df <- rbind(df_lh, df_rh)
df$hemi <- as.factor(df$hemi)
pred <- rbind(pred_lh, pred_rh)
pred$hemi <- as.factor(pred$hemi)
p1 <- ggplot(data = df, aes_string(x = smooth_var,y = thisResp, group = "hemi")) +
geom_ribbon(data = pred,aes_string(x = smooth_var , ymin = "selo",ymax = "sehi", fill = "hemi"),alpha = .5, linetype = 0) +
geom_line(data = pred,aes_string(x = smooth_var, y = "fit", color = "hemi"),size = line_size) +
scale_fill_manual(values = c("Left" = custom_color, "Right" = custom_color_rh)) +
scale_color_manual(values = c("Left" = custom_color, "Right" = custom_color_rh)) +
labs(title = plabels) +
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.title.y = element_blank()) # will add labels after plotting all models together
} else {
if (any(class(modobj)=="gam")) {
model <- modobj
} else if (class(modobj$gam)=="gam") {
model <- modobj$gam
} else {
stop("Can't find a gam object to plot")
}
s<-summary(model)
## Generate custom line plot
np <- 100 #number of predicted values
df = model$model
theseVars <- attr(model$terms,"term.labels")
varClasses <- attr(model$terms,"dataClasses")
thisResp <- as.character(model$terms[[2]])
# No interaction variable, just produce a single line plot
thisPred <- data.frame(init = rep(0,np))
for (v in c(1:length(theseVars))) {
thisVar <- theseVars[[v]]
thisClass <- varClasses[thisVar]
if (thisVar == smooth_var) {
thisPred[,smooth_var] = seq(min(df[,smooth_var],na.rm = T),max(df[,smooth_var],na.rm = T), length.out = np)
} else {
switch (thisClass,
"numeric" = {thisPred[,thisVar] = median(df[,thisVar])},
"factor" = {thisPred[,thisVar] = levels(df[,thisVar])[[1]]},
"ordered" = {thisPred[,thisVar] = levels(df[,thisVar])[[1]]}
)
}
}
pred <- thisPred %>% select(-init)
p<-data.frame(predict(model,pred,se.fit = T))
pred <- cbind(pred,p)
pred$selo <- pred$fit - 2*pred$se.fit
pred$sehi <- pred$fit + 2*pred$se.fit
pred[,group_var] = NA
pred[,thisResp] = 1
p1 <- ggplot(data = df, aes_string(x = smooth_var,y = thisResp)) +
geom_ribbon(data = pred,aes_string(x = smooth_var, ymin = "selo", ymax = "sehi"),alpha = .5, linetype = 0, fill = custom_color) +
geom_line(data = pred,aes_string(x = smooth_var, y = "fit"),size = line_size, color = custom_color) +
labs(title = plabels) +
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.title.y = element_blank()) # will add labels after plotting all models together
}
if (derivative_plot == T) {
# We will add a bar that shows where the derivative is significant.
# First make some adjustments to the line plot.
p1<- p1+theme(text = element_text(size=this_font_size),
axis.text = element_text(size = this_font_size),
axis.title.y = element_text(size = this_font_size),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
legend.text = element_text(size = this_font_size),
legend.title = element_text(size = this_font_size),
axis.title = element_text(size = this_font_size),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(fill = "transparent",colour = NA),
plot.background = element_rect(fill = "transparent",colour = NA),
plot.margin = unit(c(.2, .2, 0, .2), "cm")) #Top, left,Bottom, right
scatter <- list(p1)
# Now add the plots using the derivative plotting function
if (any(grepl(x = row.names(s$s.table),pattern = ":") & grepl(x=row.names(s$s.table),pattern = int_var))) {
# Factor levels separately if there is an interaction in the model.
f<-formula(model) # current formula
fterms <- terms(f)
fac <- attr(fterms, "factors")
idx <- which(as.logical(colSums(fac[grep(x=row.names(fac),pattern = int_var),])))
new_terms <- drop.terms(fterms, dropx = idx, keep.response = TRUE)
new_formula <- formula(new_terms) # Formula without any interaction terms in the model.
#add derivative gradients for each level of the factor
num_levels <- length(levels(df[,int_var]))
level_colors <- suppressWarnings(RColorBrewer::brewer.pal(num_levels,"Set1")) #use the same palette as the line plot
plotlist = vector(mode = "list",length = num_levels+1) # we will be building a list of plots
plotlist[1] = scatter # first the scatter plot
for (fcount in 1:num_levels) {
this_level <- levels(df[,int_var])[fcount]
df$subset <- df[,int_var] == this_level
df$group_var <- df[,group_var]
this_mod <- gam(formula = new_formula,data = df,subset = subset,random=list(group_var=~1))
# this_d <- get_derivs_and_plot(modobj = this_mod,smooth_var = smooth_var,low_color = "white",hi_color = level_colors[fcount])
this_d <- get_derivs_and_plot(modobj = this_mod,smooth_var = smooth_var,low_color = "white",hi_color = level_colors[fcount])
if (fcount != num_levels & fcount != 1){
# get rid of redundant junk
this_d$theme$axis.title = element_blank()
this_d$theme$axis.text.x = element_blank()
this_d$theme$axis.ticks=element_blank()
this_d$theme$legend.background=element_blank()
this_d$theme$legend.box.background = element_blank()
this_d$theme$legend.key = element_blank()
this_d$theme$legend.title = element_blank()
this_d$theme$legend.text = element_blank()
}
if (fcount == 1) {
this_d$theme$axis.title = element_blank()
this_d$theme$axis.text.x = element_blank()
this_d$theme$axis.ticks=element_blank()
this_d$theme$legend.background=element_blank()
this_d$theme$legend.box.background = element_blank()
this_d$theme$legend.key = element_blank()
this_d$theme$legend.text = element_blank()
}
if (fcount == num_levels) {
this_d$theme$legend.background=element_blank()
this_d$theme$legend.box.background = element_blank()
this_d$theme$legend.key = element_blank()
this_d$theme$legend.title = element_blank()
}
this_d$labels$fill=NULL
plotlist[fcount+1] = list(this_d)
}
pg<-plot_grid(rel_heights = c(16*num_levels,rep(num_levels,num_levels-1),3*num_levels),plotlist = plotlist,align = "v",axis = "lr",ncol = 1)
final_plot <- pg
print(final_plot)
} else {
# No need to split
d1 <- get_derivs_and_plot(modobj = modobj,smooth_var = smooth_var)
scatter <- list(p1)
bar <- list(d1)
allplots <- c(scatter,bar)
pg<-plot_grid(rel_heights = c(16,3),plotlist = allplots,align = "v",axis = "lr",ncol = 1)
final_plot <- pg
print(final_plot)
}
} else {
# No derivative plot
p1<- p1+theme(text = element_text(size=font_size),
axis.text = element_text(size = font_size),
legend.text = element_text(size = font_size),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
plot.background = element_blank())
final_plot <- p1
print(final_plot)
}
if (check_diagnostics == T) {
cp <- check(b,
a.qq = list(method = "tnorm",
a.cipoly = list(fill = "light blue")),
a.respoi = list(size = 0.5),
a.hist = list(bins = 10))
print(cp)
}
return(final_plot)
}
# function to implement visualize_model_hemis for all tracts
visualize_tractprofiles <- function(tract, scalar, df) {
if(str_detect(tract, "Forceps")) {
tract_data <- df %>% filter(tractID == tract)
tract_model <- gamm(as.formula(sprintf("%s ~ s(nodeID, k=24, fx = F)", scalar)), random = list(subjectID=~1), data=tract_data)
print(paste("Finished fitting GAM for", tract))
plot <- visualize_model_hemis(tract_model, both_hemis = F, smooth_var = "nodeID", int_var = NULL, group_var="subjectID", plabels = NULL,
check_diagnostics = F,derivative_plot = F, custom_color = "#831D69FF") +
labs(title = tract) + theme(plot.title = element_text(size=14, hjust=0.5))
print(paste("Finished plotting", tract))
} else {
tract_data_lh <- df %>% filter(tractID == tract) %>% filter(str_detect(tract_hemi, "Left"))
tract_data_rh <- df %>% filter(tractID == tract) %>% filter(str_detect(tract_hemi, "Right"))
tract_model_lh <- gamm(as.formula(sprintf("%s ~ s(nodeID, k=24, fx = F)", scalar)), random = list(subjectID=~1), data=tract_data_lh)
print(paste("Finished fitting GAM for left", tract))
tract_model_rh <- gamm(as.formula(sprintf("%s ~ s(nodeID, k=24, fx = F)", scalar)), random = list(subjectID=~1), data=tract_data_rh)
print(paste("Finished fitting GAM for right", tract))
plot <- visualize_model_hemis(tract_model_lh$gam, tract_model_rh$gam, both_hemis = T, smooth_var = "nodeID", int_var = NULL, group_var="subjectID", plabels = NULL,
check_diagnostics = F,derivative_plot = F, custom_color = "#831D69FF", custom_color_rh = "#FA9BAC") +
labs(title = tract) + theme(plot.title = element_text(size=14, hjust=0.5))
print(paste("Finished plotting", tract))
}
return(plot)
}
# function for plotting tractprofiles mean and sd FA/MD
plot_mean_sd <- function(tract, scalar, custom_color, custom_color_rh=NULL) {
df_profiles <- get(paste0("tract_profiles_", scalar, "_summary"))
df <- df_profiles %>% filter(tractID == tract)
if(is.null(custom_color_rh)) {
plot <- ggplot(data = df, aes_string(x = "nodeID", y = scalar)) +
geom_ribbon(data = df, aes(x = nodeID , ymin = ymin ,ymax = ymax), alpha = .5, linetype = 0, fill= custom_color) +
geom_line(data = df, aes(x = nodeID, y = mean_scalar), size = line_size, color = custom_color) +
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.title.y = element_blank())
} else {
plot <- ggplot(data = df, aes(x = nodeID, y = mean_scalar, group = hemi)) +
geom_ribbon(data = df, aes(x = nodeID , ymin = ymin ,ymax = ymax, fill = hemi), alpha = .5, linetype = 0) +
geom_line(data = df, aes(x = nodeID, y = mean_scalar, color = hemi),size = 1) +
scale_fill_manual(values = c("Left" = custom_color, "Right" = custom_color_rh)) +
scale_color_manual(values = c("Left" = custom_color, "Right" = custom_color_rh)) +
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.title.y = element_blank(),
plot.title = element_text(size=14, hjust=0.5)) + labs(title = tract)
}
return(plot)
}
arrange_by_orientation <- function(list_figures, age_effect_type) {
AP_plots <- ggarrange(list_figures$`Anterior Thalamic Radiation`,
list_figures$`Cingulum Cingulate`,
list_figures$`Inferior Fronto-occipital Fasciculus`,
list_figures$`Inferior Longitudinal Fasciculus`,
list_figures$`Superior Longitudinal Fasciculus`, ncol=5, nrow = 1)
AP_plots_final <- annotate_figure(AP_plots, top = text_grob("Anterior - Posterior",
color = "black", face = "bold", size = 14))
AP_frontotemp_plots <- ggarrange(list_figures$`Arcuate Fasciculus`,
list_figures$`Uncinate Fasciculus`, ncol=2, nrow = 1)
AP_frontotemp_plots_final <- annotate_figure(AP_frontotemp_plots, top = text_grob("Anterior - Posterior (Frontal - Temporal)",
color = "black", face = "bold", size = 14))
RL_plots <- ggarrange(list_figures$`Forceps Major`, list_figures$`Forceps Minor`,
ncol=2, nrow = 1)
RL_plots_final <- annotate_figure(RL_plots, top = text_grob("Right - Left",
color = "black", face = "bold", size = 14))
SI_plots <- ggarrange(list_figures$`Corticospinal Tract`, list_figures$`Posterior Arcuate`,
list_figures$`Vertical Occipital Fasciculus`, common.legend=TRUE, legend=c("right"),
ncol=3, nrow = 1)
SI_plots_final <- annotate_figure(SI_plots, top = text_grob("Superior - Inferior",
color = "black", face = "bold", size = 14))
tractprofiles_plot <- ggarrange(AP_plots_final, ggarrange(AP_frontotemp_plots_final, RL_plots_final, ncol = 2), SI_plots_final, nrow = 3)
tractprofiles_plot_final <- annotate_figure(tractprofiles_plot, top = text_grob(age_effect_type,
color = "black", face = "italic", size = 18))
return(tractprofiles_plot_final)
}Fractional Anisotropy
Mean FA across subjects
tract_profiles_fa <- tract_profiles %>% select(-dti_md)
tract_profiles_fa_summary <- tract_profiles_fa %>%
group_by(tractID, tract_hemi, nodeID, hemi) %>%
summarise(mean_scalar = mean(dti_fa, na.rm=T),
sd = sd(dti_fa, na.rm=T),
ymin = mean_scalar - sd,
ymax = mean_scalar + sd)
# plot
mean_sd_fa <- lapply(unique(tract_profiles$tractID), plot_mean_sd, scalar="fa", custom_color = "#831D69FF", custom_color_rh = "#FA9BAC")
names(mean_sd_fa) <- unique(tract_profiles$tractID)
# arrange by AP, AP_frontal_temporal, RL (forceps), and SI
x.grob <- textGrob("Position on Tract (Node ID)",
gp=gpar(col="black", fontsize=14))
y.grob <- textGrob("Fractional Anisotropy",
gp=gpar(col="black", fontsize=14), rot=90)
space.grob <- textGrob(" ",
gp=gpar(col="black", fontsize=28), rot=90)
tractprofiles_fa_mean_final <- arrange_by_orientation(mean_sd_fa, "Mean +/- 1 SD")
grid.arrange(arrangeGrob(tractprofiles_fa_mean_final, left = y.grob, bottom = x.grob, right = space.grob))
#ggsave("/Users/audluo/PennLINC/luowm_local/output/tract_profiles_testing/HCPD/tract_profiles_fa_mean.png", grid.arrange(arrangeGrob(tractprofiles_fa_mean_final, left = y.grob, bottom = x.grob, right = space.grob)), width = 15, height = 10, units = "in")FA tract profiles modeled using GAMs for each tract:
gamm(FA ~ s(nodeID, k=24, fx = F), random = list(subjectID=~1), data=tract_data)
# plot by AP, AP_frontal_temporal, RL (forceps), and SI
tractprofiles_fa_final <- arrange_by_orientation(tractprofiles_fa_plots, "")
grid.arrange(arrangeGrob(tractprofiles_fa_final, left = y.grob, bottom = x.grob, right = space.grob))
Mean Diffusivity
Mean MD across subjects
tract_profiles_dti_md <- tract_profiles %>% select(-dti_fa)
sqrt_n <- length(unique(tract_profiles_dti_md$subjectID))
tract_profiles_dti_md_summary <- tract_profiles_dti_md %>%
group_by(tractID, tract_hemi, nodeID, hemi) %>%
summarise(mean_scalar = mean(dti_md, na.rm=T),
sd = sd(dti_md, na.rm=T),
#se = sd/sqrt_n,
#ymin = mean_scalar - se,
#ymax = mean_scalar + se,
ymin = mean_scalar - sd,
ymax = mean_scalar + sd)
# plot
mean_sd_dti_md <- lapply(unique(tract_profiles$tractID), plot_mean_sd, scalar="dti_md", custom_color = "#831D69FF", custom_color_rh = "#FA9BAC")
names(mean_sd_dti_md) <- unique(tract_profiles$tractID)
# arrange by AP, AP_frontal_temporal, RL (forceps), and SI
x.grob <- textGrob("Position on Tract (Node ID)",
gp=gpar(col="black", fontsize=14))
y.grob <- textGrob("Mean Diffusivity",
gp=gpar(col="black", fontsize=14), rot=90)
space.grob <- textGrob(" ",
gp=gpar(col="black", fontsize=28), rot=90)
# arrange by AP, AP_frontal_temporal, RL (forceps), and SI
tractprofiles_md_mean_final <- arrange_by_orientation(mean_sd_dti_md, "Mean +/- 1 SD")
grid.arrange(arrangeGrob(tractprofiles_md_mean_final, left = y.grob, bottom = x.grob, right = space.grob))
#ggsave("/Users/audluo/PennLINC/luowm_local/output/tract_profiles_testing/HCPD/tract_profiles_dti_md_mean.png", grid.arrange(arrangeGrob(tractprofiles_md_mean_final, left = y.grob, bottom = x.grob, right = space.grob)), width = 15, height = 10, units = "in")MD tract profiles modeled using GAMs for each tract:
gamm(MD ~ s(nodeID, k=24, fx = F), random = list(subjectID=~1), data=tract_data)
# plot by AP, AP_frontal_temporal, RL (forceps), and SI
x.grob <- textGrob("Position on Tract (Node ID)",
gp=gpar(col="black", fontsize=14))
y.grob <- textGrob("Mean Diffusivity",
gp=gpar(col="black", fontsize=14), rot=90)
space.grob <- textGrob(" ",
gp=gpar(col="black", fontsize=28), rot=90)
# plot by AP, AP_frontal_temporal, RL (forceps), and SI
tractprofiles_md_final <- arrange_by_orientation(tractprofiles_md_plots, "")
grid.arrange(arrangeGrob(tractprofiles_md_final, left = y.grob, bottom = x.grob, right = space.grob))
2. Age effects
- Age effect was computed in 3 different ways:
- Using nodewise GAMs to compute delta adjusted rsq for the age term
- Using tractwise GAMs, and estimating posterior fitted values from the simulated GAM posterior distribution for each node age each age. Then computing the percent change in FA/MD for each node for each draw of the max and min age
- Same as 2, but using normalized posterior percent change
# function for computing credible intervals for percent change
compute_credible_intervals <- function(age_effect_df) {
df <- age_effect_df %>% select(contains("draw"))
median_age_effect <- apply(df, 1, function(x){median(x)})
age_effect.credible.intervals <- apply(df, 1, function(x){quantile(x, probs = c(0.025, 0.975))}) #compute the credible interval for the percent change at each node
age_effect.credible.intervals <- as.data.frame(t(age_effect.credible.intervals))
credible.intervals <- cbind(age_effect_df$nodeID, age_effect_df$tract_hemi, age_effect_df$tractID, age_effect_df$hemi, age_effect.credible.intervals)
credible.intervals <- cbind(credible.intervals, median_age_effect)
colnames(credible.intervals) <- c("nodeID","tract_hemi", "tractID", "hemi","lower","upper", "median_percent_change")
credible.intervals$median_percent_change <- credible.intervals$median_percent_change*100
credible.intervals$lower <- credible.intervals$lower*100
credible.intervals$upper <- credible.intervals$upper*100
return(credible.intervals)
}
# function for plotting age effect
plot_age_effect <- function(tract, age_effect_df, age_effect_type, custom_color, custom_color_rh=NULL) {
df <- age_effect_df %>% filter(tractID == tract)
if(age_effect_type == "adj_rsq") { # delta adj rsq
# NA out color/fill aes if adj rsq = 0 or if p-value doesn't survive FDR correction (makes the color gray)
includes_zero <- which(df$s_age.delta.adj.rsq_signed==0 | df$s_age.p.value.fdr > 0.05)
if(str_detect(tract, "Forceps")) {
df$tract_hemi[includes_zero] <- NA
plot <- ggplot(data = df, aes(x = nodeID, y = s_age.delta.adj.rsq_signed, colour = tract_hemi, fill = tract_hemi)) +
geom_point(shape = 21, size=2, alpha = 0.6) +
scale_colour_manual(values = custom_color, na.value = "grey50") +
scale_fill_manual(values = custom_color, na.value = "grey50") +
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.title.y = element_blank(),
plot.title = element_text(size=14, hjust=0.5)) + labs(title = tract)
} else {
df$hemi[includes_zero] <- NA
plot <- ggplot(data = df, aes(x = nodeID, y = s_age.delta.adj.rsq_signed, group = hemi, colour = hemi, fill = hemi)) +
geom_point(shape = 21, size=2, alpha = 0.8) +
scale_colour_manual(values = c("Left" = custom_color, "Right" = custom_color_rh), na.value = "grey50") +
scale_fill_manual(values = c("Left" = custom_color, "Right" = custom_color_rh), na.value = "grey50") +
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.title.y = element_blank(),
plot.title = element_text(size=14, hjust=0.5)) + labs(title = tract)
}
} else if (age_effect_type == "percent_change") { # percent change or normalized percent change
# NA out hemi if median_percent_change = 0 or if lower/upper CI contains 0 (makes the color gray)
includes_zero <- which(df$lower <= 0 & df$upper >= 0)
if(str_detect(tract, "Forceps")) {
df$tractID[includes_zero] <- NA
plot <- ggplot(data = df, aes(x = nodeID, y = median_percent_change, colour = tractID, fill = tractID)) +
geom_point(shape = 21, size=2, alpha = 0.6) +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.2,
position=position_dodge(0.05), alpha = 0.8) +
scale_colour_manual(values = custom_color, na.value = "grey50") +
scale_fill_manual(values = custom_color, na.value = "grey50") +
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.title.y = element_blank(),
plot.title = element_text(size=14, hjust=0.5)) + labs(title = tract)
} else {
df$hemi[includes_zero] <- NA
plot <- ggplot(data = df, aes(x = nodeID, y = median_percent_change, group = hemi, colour = hemi, fill = hemi)) +
geom_point(shape = 21, size=2, alpha = 0.6) +
geom_errorbar(aes(ymin=lower, ymax=upper, colour = hemi), width=.2,
position=position_dodge(0.05), alpha = 0.8) +
scale_colour_manual(values = c("Left" = custom_color, "Right" = custom_color_rh), na.value = "grey50") +
scale_fill_manual(values = c("Left" = custom_color, "Right" = custom_color_rh), na.value = "grey50") +
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.title.y = element_blank(),
plot.title = element_text(size=14, hjust=0.5)) + labs(title = tract)
}
} else {
print("Provide valid age effect type")
}
return(plot)
}###############################################
### load nodewise GAM results: adjusted Rsq ###
###############################################
# load nodewise GAM results for dti_fa and dti_md
gam_age_fa <- read.csv("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_fa/GAMresults.dti_fa.age.csv")
gam_age_fa <- gam_age_fa %>% select(-X)
gam_age_md <- read.csv("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_md/GAMresults.dti_md.age.csv")
gam_age_md <- gam_age_md %>% select(-X)
# determine correct age effect sign for delta adjusted Rsq using equivalent linear model
# load linear models
lm_age_fa <- read.csv("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_fa/LinearModelresults.dti_fa.age.csv")
lm_age_md <- read.csv("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_md/LinearModelresults.dti_md.age.csv")
# take absolute value of delta adj rsq
gam_age_fa$s_age.delta.adj.rsq_signed <- abs(gam_age_fa$s_age.delta.adj.rsq)
gam_age_md$s_age.delta.adj.rsq_signed <- abs(gam_age_md$s_age.delta.adj.rsq)
# apply the sign of linear model age estimate (t-value) to GAM delta adj rsq
gam_age_fa$s_age.delta.adj.rsq_signed <- gam_age_fa$s_age.delta.adj.rsq_signed * ifelse(lm_age_fa$age.estimate >= 0, 1, -1)
gam_age_md$s_age.delta.adj.rsq_signed <- gam_age_md$s_age.delta.adj.rsq_signed * ifelse(lm_age_md$age.estimate >= 0, 1, -1)
# How many nodes have significant age effect after FDR?
dti_fa_significant_count <- length(which(gam_age_fa$s_age.p.value.fdr < 0.05))
md_significant_count <- length(which(gam_age_md$s_age.p.value.fdr < 0.05))
total_nodes = 2200
print(paste0("Nodewise GAMs - dti_fa: ", dti_fa_significant_count, " nodes ", round(dti_fa_significant_count / total_nodes * 100, 2), "% with significant age effect after FDR correction"))## [1] "Nodewise GAMs - dti_fa: 1087 nodes 49.41% with significant age effect after FDR correction"print(paste0("Nodewise GAMs - dti_md: ", md_significant_count, " nodes ", round(md_significant_count / total_nodes * 100, 2), "% with significant age effect after FDR correction"))## [1] "Nodewise GAMs - dti_md: 2191 nodes 99.59% with significant age effect after FDR correction"# Add tractID and nodeID to gam_age_fa
gam_age_fa$tractID <- all_subjects$tractID[c(1:2200)]
gam_age_fa$tract_hemi <- all_subjects$tract_hemi[c(1:2200)]
gam_age_fa$hemi <- all_subjects$hemi[c(1:2200)]
gam_age_fa$nodeID <- all_subjects$nodeID[c(1:2200)]
# Add tractID and nodeID to gam_age_md
gam_age_md$tractID <- all_subjects$tractID[c(1:2200)]
gam_age_md$tract_hemi <- all_subjects$tract_hemi[c(1:2200)]
gam_age_md$hemi <- all_subjects$hemi[c(1:2200)]
gam_age_md$nodeID <- all_subjects$nodeID[c(1:2200)]
############################################################
### load tractwise GAM results: posterior percent change ###
############################################################
percent_change_fa <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_fa/posterior_percentchange/HCPD_posterior_percentchange.RData")
percent_change_fa <- percent_change_fa %>% mutate(tractID = gsub("_"," ", tract_hemi)) %>%
mutate(hemi = ifelse(grepl("Left", tract_hemi), "Left", "Right"))
percent_change_fa$tractID <- gsub("Left |Right ", "", percent_change_fa$tractID)
percent_change_md <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_md/posterior_percentchange/HCPD_posterior_percentchange.RData")
percent_change_md <- percent_change_md %>% mutate(tractID = gsub("_"," ", tract_hemi)) %>%
mutate(hemi = ifelse(grepl("Left", tract_hemi), "Left", "Right"))
percent_change_md$tractID <- gsub("Left |Right ", "", percent_change_md$tractID)
#######################################################################
### load tractwise GAM results: normalized posterior percent change ###
#######################################################################
norm_percent_change_fa <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_fa/posterior_percentchange/HCPD_posterior_percentchange_normalized.RData")
norm_percent_change_fa <- norm_percent_change_fa %>% mutate(tractID = gsub("_"," ", tract_hemi)) %>%
mutate(hemi = ifelse(grepl("Left", tract_hemi), "Left", "Right"))
norm_percent_change_fa$tractID <- gsub("Left |Right ", "", norm_percent_change_fa$tractID)
norm_percent_change_md <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_md/posterior_percentchange/HCPD_posterior_percentchange_normalized.RData")
norm_percent_change_md <- norm_percent_change_md %>% mutate(tractID = gsub("_"," ", tract_hemi)) %>%
mutate(hemi = ifelse(grepl("Left", tract_hemi), "Left", "Right"))
norm_percent_change_md$tractID <- gsub("Left |Right ", "", norm_percent_change_md$tractID)
#########################################################
### tractwise GAM results: compute credible intervals ###
#########################################################
ci_pc_fa <- compute_credible_intervals(percent_change_fa)
ci_pc_md <- compute_credible_intervals(percent_change_md)
ci_norm_pc_fa <- compute_credible_intervals(norm_percent_change_fa)
ci_norm_pc_md <- compute_credible_intervals(norm_percent_change_md)
######################################################
### Fix tract orientations for all age effect df's ###
######################################################
fix_tract_orientations <- function(df) {
# flip tracts
to_reorient <- df %>%
group_by(tractID) %>%
mutate(nodeID = ifelse(tractID %in% c("Cingulum Cingulate", "Inferior Fronto-occipital Fasciculus", "Inferior Longitudinal Fasciculus", "Uncinate Fasciculus", "Corticospinal Tract"), max(nodeID) - nodeID, nodeID))
# label main orientation
reoriented <- to_reorient %>%
mutate(main_orientation = case_when(
tractID %in% c("Anterior Thalamic Radiation", "Cingulum Cingulate", "Inferior Fronto-occipital Fasciculus",
"Inferior Longitudinal Fasciculus", "Superior Longitudinal Fasciculus") ~ "AP",
tractID %in% c("Arcuate Fasciculus", "Uncinate Fasciculus") ~ "AP_frontal_temporal",
tractID %in% c("Forceps Minor", "Forceps Major") ~ "RL",
tractID %in% c("Corticospinal Tract", "Posterior Arcuate","Vertical Occipital Fasciculus") ~ "SI",
TRUE ~ NA_character_
))
reoriented$main_orientation <- as.factor(reoriented$main_orientation)
return(reoriented)
}
gam_age_fa <- fix_tract_orientations(gam_age_fa)
gam_age_md <- fix_tract_orientations(gam_age_md)
ci_pc_fa <- fix_tract_orientations(ci_pc_fa)
ci_pc_md <- fix_tract_orientations(ci_pc_md)
ci_norm_pc_fa <- fix_tract_orientations(ci_norm_pc_fa)
ci_norm_pc_md <- fix_tract_orientations(ci_norm_pc_md)
#write.csv(gam_age_fa, "/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_fa/GAMresults.dti_fa.age_reoriented.csv")
#write.csv(gam_age_md, "/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_md/GAMresults.dti_md.age_reoriented.csv")
#saveRDS(ci_pc_fa, "/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_fa/posterior_percentchange/HCPD_posterior_percentchange_reoriented.RData")
#saveRDS(ci_pc_md, "/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_md/posterior_percentchange/HCPD_posterior_percentchange_reoriented.RData")
#saveRDS(ci_norm_pc_fa, "/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_fa/posterior_percentchange/HCPD_posterior_percentchange_normalized_reoriented.RData")
#saveRDS(ci_norm_pc_md, "/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_md/posterior_percentchange/HCPD_posterior_percentchange_normalized_reoriented.RData")Fractional Anisotropy
1. Nodewise GAMs: Delta Adjusted Rsq
- Delta adjust Rsq was computed using nodewise gam
gam(fa_node_i_tract_i ~ s(age) + covariates))
# plot
age_rsq_fa_plots <- lapply(unique(gam_age_fa$tractID), plot_age_effect, age_effect_df = gam_age_fa, age_effect_type = "adj_rsq", custom_color = "#831D69FF", custom_color_rh = "#FA9BAC")
names(age_rsq_fa_plots) <- unique(gam_age_fa$tractID)
# arrange by AP, AP_frontal_temporal, RL (forceps), and SI
x.grob <- textGrob("Position on Tract (Node ID)",
gp=gpar(col="black", fontsize=14))
y.grob <- textGrob("Fractional Anisotropy Delta Adj Rsq",
gp=gpar(col="black", fontsize=14), rot=90)
space.grob <- textGrob(" ",
gp=gpar(col="black", fontsize=14), rot=90)
age_rsq_fa_plots_final <- arrange_by_orientation(age_rsq_fa_plots, "FA Delta Adjusted Rsq")
grid.arrange(arrangeGrob(age_rsq_fa_plots_final, left = y.grob, bottom = x.grob, right = space.grob))
2. Tractwise GAMs: Posterior percent change
- GAM used:
gam(tract_i_fa ~ s(age, k = 3, fx = T, bs = 'cr') + s(nodeID, k = 24, fx = F, bs = 'cr') + ti(age, nodeID, k=c(3, 24), fx = F) + s(subjectID, bs = 're’) , method = “REML”) - Estimated posterior fitted values from simulated GAM posterior distribution for each node at each age
- Computed percent change for each node for each draw of the max and
min age:
(Y_old - Y_young)/(Y_young) - Computed credible intervals on the node-wise percent change from the posterior percent change values
percent_change_fa_plots <- lapply(unique(all_subjects$tractID), plot_age_effect, age_effect_df = ci_pc_fa, age_effect_type = "percent_change", custom_color = "#831D69FF", custom_color_rh = "#FA9BAC")
names(percent_change_fa_plots) <- unique(all_subjects$tractID)
# arrange by AP, AP_frontal_temporal, RL (forceps), and SI
x.grob <- textGrob("Position on Tract (Node ID)",
gp=gpar(col="black", fontsize=14))
y.grob <- textGrob("Fractional Anisotropy Percent Change",
gp=gpar(col="black", fontsize=14), rot=90)
space.grob <- textGrob(" ",
gp=gpar(col="black", fontsize=14), rot=90)
age_pc_fa_plots_final <- arrange_by_orientation(percent_change_fa_plots, "FA Percent Change")
grid.arrange(arrangeGrob(age_pc_fa_plots_final, left = y.grob, bottom = x.grob, right = space.grob))
3. Tractwise GAMs: Normalized posterior percent change
- Same as
2but computed normalized posterior percent change instead:(Y_old - Y_young)/(mean_tract_i_node_i_fa)
norm_percent_change_fa_plots <- lapply(unique(all_subjects$tractID), plot_age_effect, age_effect_df = ci_norm_pc_fa, age_effect_type = "percent_change", custom_color = "#831D69FF", custom_color_rh = "#FA9BAC")
names(norm_percent_change_fa_plots) <- unique(all_subjects$tractID)
# arrange by AP, AP_frontal_temporal, RL (forceps), and SI
x.grob <- textGrob("Position on Tract (Node ID)",
gp=gpar(col="black", fontsize=14))
y.grob <- textGrob("Fractional Anisotropy Normalized Percent Change",
gp=gpar(col="black", fontsize=14), rot=90)
space.grob <- textGrob(" ",
gp=gpar(col="black", fontsize=14), rot=90)
age_norm_pc_plots_final <- arrange_by_orientation(norm_percent_change_fa_plots, "FA Normalized Percent Change")
grid.arrange(arrangeGrob(age_norm_pc_plots_final, left = y.grob, bottom = x.grob, right = space.grob))
Mean Diffusivity
1. Nodewise GAMs: Delta Adjusted Rsq
- Delta adjust Rsq was computed using nodewise gam
gam(md_node_i_tract_i ~ s(age) + covariates))
# plot
age_rsq_md_plots <- lapply(unique(gam_age_md$tractID), plot_age_effect, age_effect_df = gam_age_md, age_effect_type = "adj_rsq", custom_color = "#831D69FF", custom_color_rh = "#FA9BAC")
names(age_rsq_md_plots) <- unique(gam_age_md$tractID)
# arrange by AP, AP_frontal_temporal, RL (forceps), and SI
x.grob <- textGrob("Position on Tract (Node ID)",
gp=gpar(col="black", fontsize=14))
y.grob <- textGrob("Mean Diffusivity Delta Adj Rsq",
gp=gpar(col="black", fontsize=14), rot=90)
space.grob <- textGrob(" ",
gp=gpar(col="black", fontsize=14), rot=90)
age_rsq_md_plots_final <- arrange_by_orientation(age_rsq_md_plots, "MD Delta Adjusted Rsq")
grid.arrange(arrangeGrob(age_rsq_md_plots_final, left = y.grob, bottom = x.grob, right = space.grob))
2. Tractwise GAMs: Posterior percent change
- GAM used:
gam(tract_i_md ~ s(age, k = 3, fx = T, bs = 'cr') + s(nodeID, k = 24, fx = F, bs = 'cr') + ti(age, nodeID, k=c(3, 24), fx = F) + s(subjectID, bs = 're’) , method = “REML”) - Estimated posterior fitted values from simulated GAM posterior distribution for each node at each age
- Computed percent change for each node for each draw of the max and
min age:
(Y_old - Y_young)/(Y_young) - Computed credible intervals on the node-wise percent change from the posterior percent change values
percent_change_md_plots <- lapply(unique(all_subjects$tractID), plot_age_effect, age_effect_df = ci_pc_md, age_effect_type = "percent_change", custom_color = "#831D69FF", custom_color_rh = "#FA9BAC")
names(percent_change_md_plots) <- unique(all_subjects$tractID)
# arrange by AP, AP_frontal_temporal, RL (forceps), and SI
x.grob <- textGrob("Position on Tract (Node ID)",
gp=gpar(col="black", fontsize=14))
y.grob <- textGrob("Mean Diffusivity Percent Change",
gp=gpar(col="black", fontsize=14), rot=90)
space.grob <- textGrob(" ",
gp=gpar(col="black", fontsize=14), rot=90)
age_pc_md_plots_final <- arrange_by_orientation(percent_change_md_plots, "MD Percent Change")
grid.arrange(arrangeGrob(age_pc_md_plots_final, left = y.grob, bottom = x.grob, right = space.grob))
3. Tractwise GAMs: Normalized posterior percent change
- Same as
2but computed normalized posterior percent change instead:(Y_old - Y_young)/(mean_tract_i_node_i_md)
norm_percent_change_md_plots <- lapply(unique(all_subjects$tractID), plot_age_effect, age_effect_df = ci_norm_pc_md, age_effect_type = "percent_change", custom_color = "#831D69FF", custom_color_rh = "#FA9BAC")
names(norm_percent_change_md_plots) <- unique(all_subjects$tractID)
# arrange by AP, AP_frontal_temporal, RL (forceps), and SI
x.grob <- textGrob("Position on Tract (Node ID)",
gp=gpar(col="black", fontsize=14))
y.grob <- textGrob("Mean Diffusivity Normalized Percent Change",
gp=gpar(col="black", fontsize=14), rot=90)
space.grob <- textGrob(" ",
gp=gpar(col="black", fontsize=14), rot=90)
age_norm_pc_plots_final <- arrange_by_orientation(norm_percent_change_md_plots, "MD Normalized Percent Change")
grid.arrange(arrangeGrob(age_norm_pc_plots_final, left = y.grob, bottom = x.grob, right = space.grob))
Compare age effects to an equivalent linear model
- From Jason Yeatman:
How do the GAM results compare to a simpler model where you just fit a linear model at each node and look at the slope of development. Are the results pretty similar?
One issue that we've been having with GAMs is that small differences in model specification have dramatic impacts on the results so I always try to confirm with the simple approach.- I checked with MD only
- Output looks very similar to delta adj rsq and percent change!
# Add tractID and nodeID to lm_age_md
lm_age_md$tractID <- all_subjects$tractID[c(1:2200)]
lm_age_md$tract_hemi <- all_subjects$tract_hemi[c(1:2200)]
lm_age_md$hemi <- all_subjects$hemi[c(1:2200)]
lm_age_md$nodeID <- all_subjects$nodeID[c(1:2200)]
lm_age_md <- fix_tract_orientations(lm_age_md)
plot_age_effect_lm <- function(tract, age_effect_df, custom_color, custom_color_rh=NULL) {
df <- age_effect_df %>% filter(tractID == tract)
# NA out color/fill aes if adj rsq = 0 or if p-value doesn't survive FDR correction (makes the color gray)
includes_zero <- which(df$age.p.value.fdr > 0.05)
if(str_detect(tract, "Forceps")) {
df$tract_hemi[includes_zero] <- NA
plot <- ggplot(data = df, aes(x = nodeID, y = age.estimate, colour = tract_hemi, fill = tract_hemi)) +
geom_point(shape = 21, size=2, alpha = 0.6) +
scale_colour_manual(values= custom_color, na.value = "grey50") + # allows for dynamic updating of tractID value
scale_fill_manual(values=custom_color, na.value = "grey50") +
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.title.y = element_blank(),
plot.title = element_text(size=14, hjust=0.5)) + labs(title = tract)
} else {
df$hemi[includes_zero] <- NA
plot <- ggplot(data = df, aes(x = nodeID, y = age.estimate, group = hemi, colour = hemi, fill = hemi)) +
geom_point(shape = 21, size=2, alpha = 0.8) +
scale_colour_manual(values = c("Left" = custom_color, "Right" = custom_color_rh), na.value = "grey50") +
scale_fill_manual(values = c("Left" = custom_color, "Right" = custom_color_rh), na.value = "grey50") +
theme(legend.position = "none",
axis.title.x = element_blank(),
axis.title.y = element_blank(),
plot.title = element_text(size=14, hjust=0.5)) + labs(title = tract)
}
return(plot)
}
lm_md_plots <- lapply(unique(all_subjects$tractID), plot_age_effect_lm, age_effect_df = lm_age_md, custom_color = "#831D69FF", custom_color_rh = "#FA9BAC")
names(lm_md_plots) <- unique(all_subjects$tractID)
# arrange by AP, AP_frontal_temporal, RL (forceps), and SI
x.grob <- textGrob("Position on Tract (Node ID)",
gp=gpar(col="black", fontsize=14))
y.grob <- textGrob("Mean Diffusivity Slope of Development (from linear model)",
gp=gpar(col="black", fontsize=14), rot=90)
space.grob <- textGrob(" ",
gp=gpar(col="black", fontsize=14), rot=90)
age_lm_plots_final <- arrange_by_orientation(lm_md_plots, "MD - Slope of Age Term (from linear model)")
grid.arrange(arrangeGrob(age_lm_plots_final, left = y.grob, bottom = x.grob, right = space.grob))
3. Age effect vs. relative node position
- Are there systematic developmental changes happening along tracts?
# relative position to end of tract: computed by difference between nodeID and 0 or 99. If nodeID < 50, then get abs diff to 0. If nodeID >= 50, then get abs diff to 99
# also, create column called tract_segment to label which half of the tract the node is in (i.e. anterior vs posterior)
compute_dist_to_cortex <- function(df) {
df <- df %>% mutate(dist_to_cortex =
case_when(
nodeID < 50 & !(tractID %in% c("Anterior Thalamic Radiation", "Corticospinal Tract")) ~ abs(nodeID - 0),
nodeID >= 50 & !(tractID %in% c("Anterior Thalamic Radiation", "Corticospinal Tract")) ~ abs(nodeID - 99),
tractID %in% c("Anterior Thalamic Radiation", "Corticospinal Tract") ~ nodeID )) %>%
mutate(tract_segment =
case_when(
nodeID < 50 & tractID %in% c("Cingulum Cingulate", "Inferior Fronto-occipital Fasciculus", "Inferior Longitudinal Fasciculus", "Superior Longitudinal Fasciculus") ~ "Anterior",
nodeID >= 50 & tractID %in% c("Cingulum Cingulate", "Inferior Fronto-occipital Fasciculus", "Inferior Longitudinal Fasciculus", "Superior Longitudinal Fasciculus") ~ "Posterior",
nodeID < 50 & tractID %in% c("Arcuate Fasciculus", "Uncinate Fasciculus") ~ "Anterior_frontal",
nodeID >= 50 & tractID %in% c("Arcuate Fasciculus", "Uncinate Fasciculus") ~ "Posterior_temporal",
nodeID < 50 & tractID %in% c("Forceps Minor", "Forceps Major") ~ "Right",
nodeID >= 50 & tractID %in% c("Forceps Minor", "Forceps Major") ~ "Left",
nodeID < 50 & tractID %in% c("Posterior Arcuate", "Vertical Occipital Fasciculus") ~ "Superior",
nodeID >= 50 & tractID %in% c("Posterior Arcuate", "Vertical Occipital Fasciculus") ~ "Inferior",
nodeID < 50 & tractID %in% c("Anterior Thalamic Radiation") ~ "Anterior",
nodeID >= 50 & tractID %in% c("Anterior Thalamic Radiation") ~ "Posterior_thalamus",
nodeID < 50 & tractID %in% c("Corticospinal Tract") ~ "Superior_motor",
nodeID >= 50 & tractID %in% c("Corticospinal Tract") ~ "Inferior_brainstem",
))
return(df)
}
plot_dist_to_cortex <- function(tract, age_effect_df, age_effect_type, custom_color_seg1, custom_color_seg2) {
df <- age_effect_df %>% filter(tractID == tract) %>% arrange(tract_segment)
if(age_effect_type == "adj_rsq") { # delta adj rsq
# NA out color/fill aes if adj rsq = 0 or if p-value doesn't survive FDR correction (makes the color gray)
includes_zero <- which(df$s_age.delta.adj.rsq_signed==0 | df$s_age.p.value.fdr > 0.05)
if(str_detect(tract, "Forceps")) {
df$tract_segment[includes_zero] <- NA
seg1 <- "Left"
seg2 <- "Right"
plot <- ggplot(data = df, aes(x = dist_to_cortex, y = s_age.delta.adj.rsq_signed, colour = tract_segment, fill = tract_segment, group = tract_segment)) +
geom_point(shape = 21, size=2, alpha = 0.7) +
scale_colour_manual(values=c(custom_color_seg1, custom_color_seg2) |> `names<-`(c(seg1, seg2)), na.value = "grey50") + # allows for dynamic updating
scale_fill_manual(values=c(custom_color_seg1, custom_color_seg2) |> `names<-`(c(seg1, seg2)), na.value = "grey50") + theme_classic() +
theme(legend.position = "bottom",
axis.title.x = element_blank(),
axis.title.y = element_blank(),
plot.title = element_text(size=12, hjust=0.5)) + labs(title = tract) + labs(title = tract, colour = "") + guides(fill=FALSE)
return(plot)
} else {
df$tract_segment[includes_zero] <- NA
seg1 <- unique(df$tract_segment)[which(!is.na(unique(df$tract_segment)))][1]
seg2 <- unique(df$tract_segment)[which(!is.na(unique(df$tract_segment)))][2]
plot_lh <- df %>% filter(str_detect(df$tract_hemi, "Left")) %>%
ggplot(aes(x = dist_to_cortex, y = s_age.delta.adj.rsq_signed, group = tract_segment, colour = tract_segment, fill = tract_segment)) +
geom_point(shape = 21, size=2, alpha = 0.7) +
scale_colour_manual(values=c(custom_color_seg1, custom_color_seg2) |> `names<-`(c(seg1, seg2)), na.value = "grey50") + # allows for dynamic updating of tractID value
scale_fill_manual(values=c(custom_color_seg1, custom_color_seg2) |> `names<-`(c(seg1, seg2)), na.value = "grey50") + theme_classic() +
theme(legend.position = "bottom",
axis.title.x = element_blank(),
axis.title.y = element_blank(),
plot.title = element_text(size=12, hjust=0.5)) +
labs(title = paste("Left", tract), colour = "") + guides(fill=FALSE)
plot_rh <- df %>% filter(str_detect(df$tract_hemi, "Right")) %>%
ggplot(aes(x = dist_to_cortex, y = s_age.delta.adj.rsq_signed, group = tract_segment, colour = tract_segment, fill = tract_segment)) +
geom_point(shape = 21, size=2, alpha = 0.7) +
scale_colour_manual(values=c(custom_color_seg1, custom_color_seg2) |> `names<-`(c(seg1, seg2)), na.value = "grey50") + # allows for dynamic updating of tractID value
scale_fill_manual(values=c(custom_color_seg1, custom_color_seg2) |> `names<-`(c(seg1, seg2)), na.value = "grey50") + theme_classic() +
theme(legend.position = "bottom",
axis.title.x = element_blank(),
axis.title.y = element_blank(),
plot.title = element_text(size=12, hjust=0.5)) + labs(title = paste("Right", tract), colour = "") + guides(fill=FALSE)
return(list(plot_lh, plot_rh))
}
} else if (age_effect_type == "percent_change") { # percent change or normalized percent change
# NA out hemi if median_percent_change = 0 or if lower/upper CI contains 0 (makes the color gray)
includes_zero <- which(df$lower <= 0 & df$upper >= 0)
if(str_detect(tract, "Forceps")) {
df$tractID[includes_zero] <- NA
seg1 <- "Left"
seg2 <- "Right"
plot <- ggplot(data = df, aes(x = dist_to_cortex, y = median_percent_change, colour = tract_segment, fill = tract_segment, group = tract_segment)) +
geom_point(shape = 21, size=2, alpha = 0.7) +
geom_errorbar(aes(ymin=lower, ymax=upper), width=1,
position=position_dodge(0.05), alpha = 0.8) +
scale_colour_manual(values=c(custom_color_seg1, custom_color_seg2) |> `names<-`(c(seg1, seg2)), na.value = "grey50") + # allows for dynamic updating
scale_fill_manual(values=c(custom_color_seg1, custom_color_seg2) |> `names<-`(c(seg1, seg2)), na.value = "grey50") + theme_classic() +
theme(legend.position = "bottom",
axis.title.x = element_blank(),
axis.title.y = element_blank(),
plot.title = element_text(size=12, hjust=0.5)) + labs(title = tract, colour = "") + guides(fill=FALSE)
return(plot)
} else {
seg1 <- unique(df$tract_segment)[1]
seg2 <- unique(df$tract_segment)[2]
df$tract_segment[includes_zero] <- NA
plot_lh <- df %>% filter(str_detect(df$tract_hemi, "Left")) %>%
ggplot(aes(x = dist_to_cortex, y = median_percent_change, group = tract_segment, colour = tract_segment, fill = tract_segment)) +
geom_point(shape = 21, size=2, alpha = 0.7) +
geom_errorbar(aes(ymin=lower, ymax=upper, colour = tract_segment), width=1,
position=position_dodge(0.05), alpha = 0.8) +
scale_colour_manual(values=c(custom_color_seg1, custom_color_seg2) |> `names<-`(c(seg1, seg2)), na.value = "grey50") + # allows for dynamic updating of tractID value
scale_fill_manual(values=c(custom_color_seg1, custom_color_seg2) |> `names<-`(c(seg1, seg2)), na.value = "grey50") + theme_classic() +
theme(legend.position = "bottom",
axis.title.x = element_blank(),
axis.title.y = element_blank(),
plot.title = element_text(size=12, hjust=0.5)) + labs(title = paste("Left", tract), colour = "") + guides(fill=FALSE)
plot_rh <- df %>% filter(str_detect(df$tract_hemi, "Right")) %>%
ggplot(aes(x = dist_to_cortex, y = median_percent_change, group = tract_segment, colour = tract_segment, fill = tract_segment)) +
geom_point(shape = 21, size=2, alpha = 0.7) +
geom_errorbar(aes(ymin=lower, ymax=upper, colour = tract_segment), width=1,
position=position_dodge(0.05), alpha = 0.8) +
scale_colour_manual(values=c(custom_color_seg1, custom_color_seg2) |> `names<-`(c(seg1, seg2)), na.value = "grey50") + # allows for dynamic updating of tractID value
scale_fill_manual(values=c(custom_color_seg1, custom_color_seg2) |> `names<-`(c(seg1, seg2)), na.value = "grey50") + theme_classic() +
theme(legend.position = "bottom",
axis.title.x = element_blank(),
axis.title.y = element_blank(),
plot.title = element_text(size=12, hjust=0.5)) + labs(title = paste("Right", tract), colour = "") + guides(fill=FALSE)
return(list(plot_lh, plot_rh))
}
} else {
print("Provide valid age effect type")
}
}
arrange_dist_to_cortex <- function(list_figures, age_effect_type) {
ATR_plots <- ggarrange(list_figures$`Anterior Thalamic Radiation`[[1]], list_figures$`Anterior Thalamic Radiation`[[2]], nrow = 2, common.legend=TRUE, legend = "bottom")
AP_plots <- ggarrange(list_figures$`Cingulum Cingulate`[[1]],
list_figures$`Inferior Fronto-occipital Fasciculus`[[1]],
list_figures$`Inferior Longitudinal Fasciculus`[[1]],
list_figures$`Superior Longitudinal Fasciculus`[[1]],
list_figures$`Cingulum Cingulate`[[2]],
list_figures$`Inferior Fronto-occipital Fasciculus`[[2]],
list_figures$`Inferior Longitudinal Fasciculus`[[2]],
list_figures$`Superior Longitudinal Fasciculus`[[2]], ncol=4, nrow = 2, common.legend=TRUE, legend = "bottom")
AP_ATR_plots <- ggarrange(ATR_plots, AP_plots, ncol=2, widths = c(1,4))
AP_plots_final <- annotate_figure(AP_ATR_plots, top = text_grob("Anterior - Posterior",
color = "black", face = "bold", size = 12))
AP_frontotemp_plots <- ggarrange(list_figures$`Arcuate Fasciculus`[[1]], list_figures$`Uncinate Fasciculus`[[1]],
list_figures$`Arcuate Fasciculus`[[2]], list_figures$`Uncinate Fasciculus`[[2]], ncol=2, nrow = 2, common.legend=TRUE, legend = "bottom")
AP_frontotemp_plots_final <- annotate_figure(AP_frontotemp_plots, top = text_grob("Anterior - Posterior (Frontal - Temporal)",
color = "black", face = "bold", size = 12))
RL_plots <- ggarrange(list_figures$`Forceps Major`, list_figures$`Forceps Minor`,
ncol=2, nrow = 1, common.legend=TRUE, legend = "bottom")
RL_plots_final <- annotate_figure(RL_plots, top = text_grob("Right - Left",
color = "black", face = "bold", size = 12))
CST_plots <- ggarrange(list_figures$`Corticospinal Tract`[[1]],list_figures$`Corticospinal Tract`[[2]], common.legend=TRUE, legend=c("bottom"), ncol=1, nrow = 2)
SI_plots <- ggarrange(list_figures$`Posterior Arcuate`[[1]],
list_figures$`Vertical Occipital Fasciculus`[[1]],
list_figures$`Posterior Arcuate`[[2]],
list_figures$`Vertical Occipital Fasciculus`[[2]], common.legend=TRUE, legend=c("bottom"), ncol=2, nrow = 2)
SI_CST_plots <- ggarrange(CST_plots, SI_plots, widths = c(1, 2))
SI_plots_final <- annotate_figure(SI_CST_plots, top = text_grob("Superior - Inferior",
color = "black", face = "bold", size = 12))
tractprofiles_plot <- ggarrange(AP_plots_final, ggarrange(AP_frontotemp_plots_final, RL_plots_final, ncol = 2), SI_plots_final, nrow = 3)
tractprofiles_plot_final <- annotate_figure(tractprofiles_plot, top = text_grob(age_effect_type,
color = "black", face = "italic", size = 18))
return(tractprofiles_plot_final)
}gam_age_fa <- compute_dist_to_cortex(gam_age_fa)
gam_age_md <- compute_dist_to_cortex(gam_age_md)
ci_pc_fa <- compute_dist_to_cortex(ci_pc_fa)
ci_pc_md <- compute_dist_to_cortex(ci_pc_md)
ci_norm_pc_fa <- compute_dist_to_cortex(ci_norm_pc_fa)
ci_norm_pc_md <- compute_dist_to_cortex(ci_norm_pc_md) custom_color_seg1 = "#4CC3FFFF"
custom_color_seg2 = "#FF9932FF"
rsq_fa_dist_plots <- lapply(unique(gam_age_fa$tractID), plot_dist_to_cortex, age_effect_df = gam_age_fa, age_effect_type = "adj_rsq", custom_color_seg1, custom_color_seg2)
names(rsq_fa_dist_plots) <- unique(gam_age_fa$tractID)
# FA delta adj rsq
x.grob <- textGrob("Distance to Cortex",
gp=gpar(col="black", fontsize=12))
y.grob <- textGrob("Fractional Anisotropy Delta Adjusted Rsq",
gp=gpar(col="black", fontsize=12), rot=90)
space.grob <- textGrob(" ",
gp=gpar(col="black", fontsize=28), rot=90)
grid.arrange(arrangeGrob(arrange_dist_to_cortex(rsq_fa_dist_plots, "FA Delta Adjusted Rsq"), left = y.grob, bottom = x.grob, right = space.grob))
#ggsave("/Users/audluo/PennLINC/luowm_local/output/tract_profiles_testing/HCPD/dist_to_cort_rsq_fa.png", grid.arrange(arrangeGrob(arrange_dist_to_cortex(rsq_fa_dist_plots, "FA Delta Adjusted Rsq"), left = y.grob, bottom = x.grob, right = space.grob)), width = 16, height = 20, units = "in")pc_fa_dist_plots <- lapply(unique(ci_pc_fa$tractID), plot_dist_to_cortex, age_effect_df = ci_pc_fa, age_effect_type = "percent_change", custom_color_seg1, custom_color_seg2)
names(pc_fa_dist_plots) <- unique(ci_pc_fa$tractID)
# FA Percent Change
y.grob <- textGrob("Fractional Anisotropy Percent Change",
gp=gpar(col="black", fontsize=12), rot=90)
arrange_dist_to_cortex(pc_fa_dist_plots, "FA Percent Change")
#ggsave("/Users/audluo/PennLINC/luowm_local/output/tract_profiles_testing/HCPD/dist_to_cort_pc_fa.png", grid.arrange(arrangeGrob(arrange_dist_to_cortex(pc_fa_dist_plots, "FA Percent Change"), left = y.grob, bottom = x.grob, right = space.grob)), width = 16, height = 20, units = "in")norm_pc_fa_dist_plots <- lapply(unique(ci_norm_pc_fa$tractID), plot_dist_to_cortex, age_effect_df = ci_norm_pc_fa, age_effect_type = "percent_change", custom_color_seg1, custom_color_seg2)
names(norm_pc_fa_dist_plots) <- unique(ci_norm_pc_fa$tractID)
# FA Normalized Percent Change
y.grob <- textGrob("Fractional Anisotropy Normalized Percent Change",
gp=gpar(col="black", fontsize=12), rot=90)
arrange_dist_to_cortex(norm_pc_fa_dist_plots, "FA Normalized Percent Change")
#ggsave("/Users/audluo/PennLINC/luowm_local/output/tract_profiles_testing/HCPD/dist_to_cort_norm_pc_fa.png", grid.arrange(arrangeGrob(arrange_dist_to_cortex(norm_pc_fa_dist_plots, "FA Normalized Percent Change"), left = y.grob, bottom = x.grob, right = space.grob)), width = 16, height = 20, units = "in")rsq_md_dist_plots <- lapply(unique(gam_age_md$tractID), plot_dist_to_cortex, age_effect_df = gam_age_md, age_effect_type = "adj_rsq", custom_color_seg1, custom_color_seg2)
names(rsq_md_dist_plots) <- unique(gam_age_md$tractID)
# MD delta adj rsq
y.grob <- textGrob("Mean Diffusivity Delta Adjusted Rsq",
gp=gpar(col="black", fontsize=12), rot=90)
arrange_dist_to_cortex(rsq_md_dist_plots, "MD Delta Adjusted Rsq")
#ggsave("/Users/audluo/PennLINC/luowm_local/output/tract_profiles_testing/HCPD/dist_to_cort_rsq_md.png", grid.arrange(arrangeGrob(arrange_dist_to_cortex(rsq_md_dist_plots, "MD Delta Adjusted Rsq"), left = y.grob, bottom = x.grob, right = space.grob)), width = 16, height = 20, units = "in")pc_md_dist_plots <- lapply(unique(ci_pc_md$tractID), plot_dist_to_cortex, age_effect_df = ci_pc_md, age_effect_type = "percent_change", custom_color_seg1, custom_color_seg2)
names(pc_md_dist_plots) <- unique(ci_pc_md$tractID)
# MD Percent Change
y.grob <- textGrob("Mean Diffusivity Percent Change",
gp=gpar(col="black", fontsize=12), rot=90)
arrange_dist_to_cortex(pc_md_dist_plots, "MD Percent Change")
#ggsave("/Users/audluo/PennLINC/luowm_local/output/tract_profiles_testing/HCPD/dist_to_cort_pc_md.png", grid.arrange(arrangeGrob(arrange_dist_to_cortex(pc_md_dist_plots, "MD Percent Change"), left = y.grob, bottom = x.grob, right = space.grob)), width = 16, height = 20, units = "in")norm_pc_md_dist_plots <- lapply(unique(ci_norm_pc_md$tractID), plot_dist_to_cortex, age_effect_df = ci_norm_pc_md, age_effect_type = "percent_change", custom_color_seg1, custom_color_seg2)
names(norm_pc_md_dist_plots) <- unique(ci_norm_pc_md$tractID)
# MD Normalized Percent Change
y.grob <- textGrob("Mean Diffusivity Normalized Percent Change",
gp=gpar(col="black", fontsize=12), rot=90)
arrange_dist_to_cortex(norm_pc_md_dist_plots, "MD Normalized Percent Change")
#ggsave("/Users/audluo/PennLINC/luowm_local/output/tract_profiles_testing/HCPD/dist_to_cort_norm_pc_md.png", grid.arrange(arrangeGrob(arrange_dist_to_cortex(norm_pc_md_dist_plots, "MD Normalized Percent Change"), left = y.grob, bottom = x.grob, right = space.grob)), width = 16, height = 20, units = "in")Correlations between MD age effect and relative position to tract end (Spearman’s, since relative position is directly related to nodeID)
# make a dataframe that has column tract, correlation, and p-value
correlate_age_effect <- function(tract, age_effect_df, age_effect_type, hemi_avg = FALSE) {
df <- age_effect_df %>% filter(tractID == tract)
if (age_effect_type == "adj_rsq") {
age_effect_var <- "s_age.delta.adj.rsq_signed"
} else {
age_effect_var <- "median_percent_change"
}
if (hemi_avg == TRUE & !str_detect(tract, "Forceps")) {
df_avg <- df %>% group_by(nodeID) %>% summarise(mean=mean(get(age_effect_var)))
names(df_avg) <- c("nodeID", age_effect_var)
df <- right_join(df[c(1:100), c("tractID", "nodeID", "main_orientation", "dist_to_cortex", "tract_segment")], df_avg, by = "nodeID")
}
if(age_effect_type == "adj_rsq") {
cor_dist_cortex <- cor.test(df$dist_to_cortex, df$s_age.delta.adj.rsq_signed, method = "spearman")$estimate
pval_dist_cortex <- cor.test(df$dist_to_cortex, df$s_age.delta.adj.rsq_signed, method = "spearman")$p.value
cor_dist_along_tract <- cor.test(df$nodeID, df$s_age.delta.adj.rsq_signed, method = "spearman")$estimate
pval_dist_along_tract <- cor.test(df$nodeID, df$s_age.delta.adj.rsq_signed, method = "spearman")$p.value
} else if(age_effect_type == "percent_change") {
cor_dist_cortex <- cor.test(df$dist_to_cortex, df$median_percent_change, method = "spearman")$estimate
pval_dist_cortex <- cor.test(df$dist_to_cortex, df$median_percent_change, method = "spearman")$p.value
cor_dist_along_tract <- cor.test(df$nodeID, df$median_percent_change, method = "spearman")$estimate
pval_dist_along_tract <- cor.test(df$nodeID, df$median_percent_change, method = "spearman")$p.value
} else {
print("Provide valid age effect type")
}
cor_df <- data.frame(tract, cor_dist_cortex, pval_dist_cortex, cor_dist_along_tract, pval_dist_along_tract)
rownames(cor_df) <- NULL
cor_df <- cor_df %>% mutate(pval_dist_cortex_sig = ifelse(pval_dist_cortex < 0.05, 1, 0)) %>%
mutate(pval_dist_along_tract_sig = ifelse(pval_dist_along_tract < 0.05, 1, 0))
cor_df <- cor_df %>% select(tract, cor_dist_cortex, pval_dist_cortex, pval_dist_cortex_sig, cor_dist_along_tract, pval_dist_along_tract, pval_dist_along_tract_sig)
return(cor_df)
}
dist_age_effect_fdr <- function(cor_df) {
cor_df$pval_dist_cortex_fdr <- p.adjust(cor_df$pval_dist_cortex, method=c("fdr"))
cor_df$pval_dist_along_tract_fdr <- p.adjust(cor_df$pval_dist_along_tract, method=c("fdr"))
cor_df <- cor_df %>% mutate(pval_dist_cortex_sig = ifelse(pval_dist_cortex_fdr < 0.05, 1, 0)) %>%
mutate(pval_dist_along_tract_sig = ifelse(pval_dist_along_tract_fdr < 0.05, 1, 0))
cor_df <- cor_df %>% select(tract, cor_dist_cortex, pval_dist_cortex, pval_dist_cortex_fdr, pval_dist_cortex_sig, cor_dist_along_tract, pval_dist_along_tract, pval_dist_along_tract_fdr, pval_dist_along_tract_sig)
return(cor_df)
}
wrapper_correlate_age_effect <- function(age_effect_df, age_effect_type, tracts_exclude, hemi_avg=FALSE) {
cor_dist <- lapply(unique(age_effect_df$tractID), correlate_age_effect, age_effect_df = age_effect_df, age_effect_type = age_effect_type, hemi_avg)
cor_dist <- bind_rows(cor_dist)
cor_dist <- dist_age_effect_fdr(cor_dist)
mean_cor <- cor_dist %>% filter(!tract %in% tracts_exclude) %>%
summarise(mean_cor_dist_cortex = mean(cor_dist_cortex),
mean_cor_dist_along_tract = mean(cor_dist_along_tract))
return(list(cor_dist, mean_cor))
}
cor_dist_adj_rsq_fa <- wrapper_correlate_age_effect(gam_age_fa, age_effect_type = "adj_rsq", tracts_exclude = c("Anterior Thalamic Radiation", "Corticospinal Tract"))
cor_dist_adj_rsq_md <- wrapper_correlate_age_effect(gam_age_md, age_effect_type = "adj_rsq", tracts_exclude = c("Anterior Thalamic Radiation", "Corticospinal Tract"))
cor_dist_pc_fa <- wrapper_correlate_age_effect(ci_pc_fa, age_effect_type = "percent_change", tracts_exclude = c("Anterior Thalamic Radiation", "Corticospinal Tract"))
cor_dist_pc_md <- wrapper_correlate_age_effect(ci_pc_md, age_effect_type = "percent_change", tracts_exclude = c("Anterior Thalamic Radiation", "Corticospinal Tract"))
cor_dist_norm_pc_fa <- wrapper_correlate_age_effect(ci_norm_pc_fa, age_effect_type = "percent_change", tracts_exclude = c("Anterior Thalamic Radiation", "Corticospinal Tract"))
cor_dist_norm_pc_md <- wrapper_correlate_age_effect(ci_norm_pc_md, age_effect_type = "percent_change", tracts_exclude = c("Anterior Thalamic Radiation", "Corticospinal Tract"))
kable(cor_dist_adj_rsq_md[[1]][,c(1:3)], caption="Correlation between Age Effect (Delta Adjusted Rsq) and Position Along Tract for MD", align="ccc") %>%
kable_styling(font_size = 20, bootstrap_options = "hover")| tract | cor_dist_cortex | pval_dist_cortex |
|---|---|---|
| Anterior Thalamic Radiation | 0.0530013 | 0.4560435 |
| Cingulum Cingulate | 0.6066929 | 0.0000000 |
| Corticospinal Tract | 0.7199970 | 0.0000000 |
| Inferior Fronto-occipital Fasciculus | 0.3273996 | 0.0000022 |
| Inferior Longitudinal Fasciculus | 0.5947144 | 0.0000000 |
| Superior Longitudinal Fasciculus | 0.8436413 | 0.0000000 |
| Arcuate Fasciculus | 0.9032580 | 0.0000000 |
| Uncinate Fasciculus | 0.2254139 | 0.0013308 |
| Forceps Minor | 0.8012963 | 0.0000000 |
| Forceps Major | 0.2788537 | 0.0049629 |
| Posterior Arcuate | 0.7500994 | 0.0000000 |
| Vertical Occipital Fasciculus | 0.9354668 | 0.0000000 |
Mean correlation across cortico-cortical tracts with all nodes and after excluding nodes from each end:
# make a dataframe that has column tract, correlation, and p-value
correlate_age_effect_exclude <- function(tract, age_effect_df, age_effect_type, nodes_exclude) {
df <- age_effect_df %>% filter(tractID == tract)
to_exclude_lowest <- seq(nodes_exclude-nodes_exclude, nodes_exclude-1, 1)
to_exclude_highest <- seq(100-nodes_exclude, 100-nodes_exclude+nodes_exclude-1, 1)
df <- df %>% filter(!(nodeID %in% to_exclude_lowest | nodeID %in% to_exclude_highest))
if(age_effect_type == "adj_rsq") {
cor_dist_cortex <- cor.test(df$dist_to_cortex, df$s_age.delta.adj.rsq_signed, method = "spearman")$estimate
pval_dist_cortex <- cor.test(df$dist_to_cortex, df$s_age.delta.adj.rsq_signed, method = "spearman")$p.value
cor_dist_along_tract <- cor.test(df$nodeID, df$s_age.delta.adj.rsq_signed, method = "spearman")$estimate
pval_dist_along_tract <- cor.test(df$nodeID, df$s_age.delta.adj.rsq_signed, method = "spearman")$p.value
} else if(age_effect_type == "percent_change") {
cor_dist_cortex <- cor.test(df$dist_to_cortex, df$median_percent_change, method = "spearman")$estimate
pval_dist_cortex <- cor.test(df$dist_to_cortex, df$median_percent_change, method = "spearman")$p.value
cor_dist_along_tract <- cor.test(df$nodeID, df$median_percent_change, method = "spearman")$estimate
pval_dist_along_tract <- cor.test(df$nodeID, df$median_percent_change, method = "spearman")$p.value
} else {
print("Provide valid age effect type")
}
cor_df <- data.frame(tract, cor_dist_cortex, pval_dist_cortex, cor_dist_along_tract, pval_dist_along_tract)
rownames(cor_df) <- NULL
cor_df <- cor_df %>% mutate(pval_dist_cortex_sig = ifelse(pval_dist_cortex < 0.05, 1, 0)) %>%
mutate(pval_dist_along_tract_sig = ifelse(pval_dist_along_tract < 0.05, 1, 0))
cor_df <- cor_df %>% select(tract, cor_dist_cortex, pval_dist_cortex, pval_dist_cortex_sig, cor_dist_along_tract, pval_dist_along_tract, pval_dist_along_tract_sig)
return(cor_df)
}
wrapper_correlate_age_effect_excl <- function(age_effect_df, age_effect_type, tracts_exclude, nodes_exclude) {
cor_dist <- lapply(unique(age_effect_df$tractID), correlate_age_effect_exclude, age_effect_df = age_effect_df, age_effect_type = age_effect_type, nodes_exclude = nodes_exclude)
cor_dist <- bind_rows(cor_dist)
cor_dist <- dist_age_effect_fdr(cor_dist)
mean_cor <- cor_dist %>% filter(!tract %in% tracts_exclude) %>%
summarise(mean_cor_dist_cortex = mean(cor_dist_cortex),
mean_cor_dist_along_tract = mean(cor_dist_along_tract))
return(list(cor_dist, mean_cor))
}
cor_dist_adj_rsq_excl_fa <- wrapper_correlate_age_effect_excl(gam_age_fa, age_effect_type = "adj_rsq", tracts_exclude = c("Anterior Thalamic Radiation", "Corticospinal Tract"), nodes_exclude = 10)
cor_dist_adj_rsq_excl_md <- wrapper_correlate_age_effect_excl(gam_age_md, age_effect_type = "adj_rsq", tracts_exclude = c("Anterior Thalamic Radiation", "Corticospinal Tract"), nodes_exclude = 10)
cor_dist_pc_excl_fa <- wrapper_correlate_age_effect_excl(ci_pc_fa, age_effect_type = "percent_change", tracts_exclude = c("Anterior Thalamic Radiation", "Corticospinal Tract"), nodes_exclude = 10)
cor_dist_pc_excl_md <- wrapper_correlate_age_effect_excl(ci_pc_md, age_effect_type = "percent_change", tracts_exclude = c("Anterior Thalamic Radiation", "Corticospinal Tract"), nodes_exclude = 10)
cor_dist_norm_pc_excl_fa <- wrapper_correlate_age_effect_excl(ci_norm_pc_fa, age_effect_type = "percent_change", tracts_exclude = c("Anterior Thalamic Radiation", "Corticospinal Tract"), nodes_exclude = 10)
cor_dist_norm_pc_excl_md <- wrapper_correlate_age_effect_excl(ci_norm_pc_md, age_effect_type = "percent_change", tracts_exclude = c("Anterior Thalamic Radiation", "Corticospinal Tract"), nodes_exclude = 10)
cor_dist_adj_rsq_excl20_fa <- wrapper_correlate_age_effect_excl(gam_age_fa, age_effect_type = "adj_rsq", tracts_exclude = c("Anterior Thalamic Radiation", "Corticospinal Tract"), nodes_exclude = 20)
cor_dist_adj_rsq_excl20_md <- wrapper_correlate_age_effect_excl(gam_age_md, age_effect_type = "adj_rsq", tracts_exclude = c("Anterior Thalamic Radiation", "Corticospinal Tract"), nodes_exclude = 20)
cor_dist_pc_excl20_fa <- wrapper_correlate_age_effect_excl(ci_pc_fa, age_effect_type = "percent_change", tracts_exclude = c("Anterior Thalamic Radiation", "Corticospinal Tract"), nodes_exclude = 20)
cor_dist_pc_excl20_md <- wrapper_correlate_age_effect_excl(ci_pc_md, age_effect_type = "percent_change", tracts_exclude = c("Anterior Thalamic Radiation", "Corticospinal Tract"), nodes_exclude = 20)
cor_dist_norm_pc_excl20_fa <- wrapper_correlate_age_effect_excl(ci_norm_pc_fa, age_effect_type = "percent_change", tracts_exclude = c("Anterior Thalamic Radiation", "Corticospinal Tract"), nodes_exclude = 20)
cor_dist_norm_pc_excl20_md <- wrapper_correlate_age_effect_excl(ci_norm_pc_md, age_effect_type = "percent_change", tracts_exclude = c("Anterior Thalamic Radiation", "Corticospinal Tract"), nodes_exclude = 20)
# make a table to compare correlations
mean_cor_compare <- data.frame(rbind(c(cor_dist_adj_rsq_md[[2]][1,1], cor_dist_adj_rsq_excl_md[[2]][1,1], cor_dist_adj_rsq_excl20_md[[2]][1,1]), c(cor_dist_pc_md[[2]][1,1], cor_dist_pc_excl_md[[2]][1,1], cor_dist_pc_excl20_md[[2]][1,1]), c(cor_dist_norm_pc_md[[2]][1,1], cor_dist_norm_pc_excl_md[[2]][1,1], cor_dist_norm_pc_excl20_md[[2]][1,1]))) %>% setNames(c("All Nodes", "Exclude 10 nodes from each end", "Exclude 20 nodes from each end"))
rownames(mean_cor_compare) <- c("Delta Adj Rsq", "Posterior Percent Change", "Normalized Posterior Percent Change")
mean_cor_compare$notes <- c("CC, IFO, ILF, UNC, forceps major drop in correlation (IFO flips direction); but SLF, ARC, forceps minor, pARC and VOF are the same or better", "forceps maj reverses direction, IFO decreases; but forceps minor, CC, ARC, ILF, pARC, SLF, UNC, VOF increase in corr or stays approx the same ", "")
kable(mean_cor_compare, caption="Mean Correlation of Age Effect and Distance to Cortex") %>%
kable_styling(font_size = 20)| All Nodes | Exclude 10 nodes from each end | Exclude 20 nodes from each end | notes | |
|---|---|---|---|---|
| Delta Adj Rsq | 0.6266836 | 0.4572489 | 0.2660743 | CC, IFO, ILF, UNC, forceps major drop in correlation (IFO flips direction); but SLF, ARC, forceps minor, pARC and VOF are the same or better |
| Posterior Percent Change | 0.5960467 | 0.5925964 | 0.4889670 | forceps maj reverses direction, IFO decreases; but forceps minor, CC, ARC, ILF, pARC, SLF, UNC, VOF increase in corr or stays approx the same |
| Normalized Posterior Percent Change | 0.6005321 | 0.5959221 | 0.4900591 |
Mean correlations that were computed after averaging age effects of left and right tracts:
cor_dist_adj_rsq_fa_hemi <- wrapper_correlate_age_effect(gam_age_fa, age_effect_type = "adj_rsq", tracts_exclude = c("Anterior Thalamic Radiation", "Corticospinal Tract"), hemi_avg = TRUE)
cor_dist_adj_rsq_md_hemi <- wrapper_correlate_age_effect(gam_age_md, age_effect_type = "adj_rsq", tracts_exclude = c("Anterior Thalamic Radiation", "Corticospinal Tract"), hemi_avg = TRUE)
cor_dist_pc_fa_hemi <- wrapper_correlate_age_effect(ci_pc_fa, age_effect_type = "percent_change", tracts_exclude = c("Anterior Thalamic Radiation", "Corticospinal Tract"), hemi_avg = TRUE)
cor_dist_pc_md_hemi <- wrapper_correlate_age_effect(ci_pc_md, age_effect_type = "percent_change", tracts_exclude = c("Anterior Thalamic Radiation", "Corticospinal Tract"), hemi_avg = TRUE)
cor_dist_norm_pc_fa_hemi <- wrapper_correlate_age_effect(ci_norm_pc_fa, age_effect_type = "percent_change", tracts_exclude = c("Anterior Thalamic Radiation", "Corticospinal Tract"), hemi_avg = TRUE)
cor_dist_norm_pc_md_hemi <- wrapper_correlate_age_effect(ci_norm_pc_md, age_effect_type = "percent_change", tracts_exclude = c("Anterior Thalamic Radiation", "Corticospinal Tract"), hemi_avg = TRUE)
# make a table to compare correlations
hemi_cor_compare <- data.frame(rbind(c(cor_dist_adj_rsq_md[[2]][1,1], cor_dist_adj_rsq_md_hemi[[2]][1,1]), c(cor_dist_pc_md[[2]][1,1], cor_dist_pc_md_hemi[[2]][1,1]), c(cor_dist_norm_pc_md[[2]][1,1], cor_dist_norm_pc_md_hemi[[2]][1,1]))) %>% setNames(c("LH/RH age effects as separate datapoints", "Averaged LH/RH age effects"))
rownames(hemi_cor_compare) <- c("Delta Adj Rsq", "Posterior Percent Change", "Normalized Posterior Percent Change")
kable(hemi_cor_compare, caption="Mean Correlation of Age Effect and Relative Positive to Tract End - Average Both Hemispheres", align = "cc") %>%
kable_styling(font_size = 20, bootstrap_options = "hover")| LH/RH age effects as separate datapoints | Averaged LH/RH age effects | |
|---|---|---|
| Delta Adj Rsq | 0.6266836 | 0.6360798 |
| Posterior Percent Change | 0.5960467 | 0.5222018 |
| Normalized Posterior Percent Change | 0.6005321 | 0.5278102 |
4. Developmental trajectories along each tract
- Plot trajectories for each node and color by nodeID. Each line corresponds to a unique node.
- Note that the labels for tract orientations correspond to
low_nodeID - high_nodeID. i.e. Anterior-posterior corresponds to anterior regions having low node IDs (blue) and posterior regions having high nodeIDs (blue) - I used nodewise GAMs and estimated zero-centered smooths using
smooth_estimates(I don’t know how to usesmooth_estimatesfor tractwise GAMs - can’t figure out how to navigate multiple smooths and a tensor product lol)
cohortfile <- cohortfile %>% select(subjectID, age, sex, mean_fd)
gam_df <- merge(tract_profiles, cohortfile, by="subjectID")
gam_df <- gam_df %>% mutate(tract_node = paste0(tract_hemi, "_", nodeID))
gam_df$sex <- as.factor(gam_df$sex)
# create df with column tract_node and another column (dti_metric) that contains the values (nrow = N * n_tract_nodes)
# fa
nodewise_gam_fa <- gam_df %>% select(dti_fa, age, sex, mean_fd, tract_node)
nodewise_smoothfits_fa <- nodewise_gam_fa %>% group_by(tract_node) %>%
do(fit = smooth_estimates(gam(dti_fa ~ s(age,k=3,fx=F),data=.))) %>% unnest(fit)
nodewise_smoothfits_fa <- nodewise_smoothfits_fa %>% mutate(tractID = gsub("Left |Right ", "", gsub("_", " ", gsub("_[0-9]+", "", tract_node)))) %>%
mutate(nodeID = str_extract(tract_node, "[0-9]+")) %>%
mutate(hemi = str_extract(tract_node, "Left|Right"))
orientations <- tract_profiles[,c("tractID", "main_orientation")]
orientations <- orientations[!duplicated(orientations),]
nodewise_smoothfits_fa <- merge(nodewise_smoothfits_fa, orientations, by = "tractID")
# md
nodewise_gam_md <- gam_df %>% select(dti_md, age, sex, mean_fd, tract_node)
nodewise_smoothfits_md <- nodewise_gam_md %>% group_by(tract_node) %>%
do(fit = smooth_estimates(gam(dti_md ~ s(age,k=3,fx=F),data=.))) %>% unnest(fit)
nodewise_smoothfits_md <- nodewise_smoothfits_md %>% mutate(tractID = gsub("Left |Right ", "", gsub("_", " ", gsub("_[0-9]+", "", tract_node)))) %>%
mutate(nodeID = str_extract(tract_node, "[0-9]+")) %>%
mutate(hemi = str_extract(tract_node, "Left|Right"))
nodewise_smoothfits_md <- merge(nodewise_smoothfits_md, orientations, by = "tractID")
# plot trajectory for each tract and color by node
# maybe compute average trajectory for each tract end?plot_tract_ends_zeroed <- function(orientation, df, scalar, mycolors, ylim1, ylim2) {
if (orientation == "CP") {
title1 <- "Peripheral"
title2 <- "Central"
plot1 <- df %>% filter(tractID != "Anterior Thalamic Radiation") %>%
filter(tractID != "Corticospinal Tract") %>%
filter(nodeID < 10 | nodeID > 90) %>% group_by(tract_node) %>%
ggplot(aes(x = age, y = est, group = tract_node)) +
geom_line(aes(colour = tractID), size = 0.8, alpha = 0.7) +
theme_classic() +
theme(text=element_text(size=12),
legend.position = "right",
plot.title = element_text(hjust=0.5)) +
scale_colour_manual(values = mycolors) +
labs(y = paste("Fitted", scalar), x = "Age", title=title1) + ylim(ylim1, ylim2)
plot2 <- df %>% filter(tractID != "Anterior Thalamic Radiation") %>%
filter(tractID != "Corticospinal Tract") %>%
filter(nodeID > 39 & nodeID < 60) %>% group_by(tract_node) %>%
ggplot(aes(x = age, y = est, group = tract_node)) +
geom_line(aes(colour = tractID), size = 0.8, alpha = 0.7) +
theme_classic() +
theme(text=element_text(size=12),
legend.position = "right",
plot.title = element_text(hjust=0.5)) +
scale_colour_manual(values = mycolors) +
labs(y = paste("Fitted", scalar), x = "Age", title=title2) + ylim(ylim1, ylim2)
} else {
if (orientation == "AP") {
title1 <- "Anterior"
title2 <- "Posterior"
} else if (orientation == "AP_frontal_temporal") {
title1 <- "Frontal"
title2 <- "Temporal"
} else if (orientation == "SI") {
title1 <- "Superior"
title2 <- "Inferior"
} else {
print("Provide valid orientation (AP, AP_frontal_temporal, SI, CP)")
}
plot1 <- df %>% filter(main_orientation == orientation) %>%
filter(tractID != "Anterior Thalamic Radiation") %>%
filter(tractID != "Corticospinal Tract") %>%
filter(nodeID < 10) %>% group_by(tract_node) %>%
ggplot( aes(x = age, y = est, group = tract_node)) +
geom_line(aes(colour = tractID), size = 0.8, alpha = 0.7) + theme_classic() +
theme(text=element_text(size=12),
legend.position = "right",
plot.title = element_text(hjust=0.5)) +
scale_colour_manual(values = mycolors) +
labs(y = paste("Fitted", scalar), x = "Age", title=title1) + ylim(ylim1, ylim2)
plot2 <- df %>% filter(main_orientation == orientation) %>%
filter(tractID != "Anterior Thalamic Radiation") %>%
filter(tractID != "Corticospinal Tract") %>%
filter(nodeID > 89) %>% group_by(tract_node) %>%
ggplot( aes(x = age, y = est, group = tract_node)) +
geom_line(aes(colour = tractID), size = 0.8, alpha = 0.7) + theme_classic() +
theme(text=element_text(size=12),
legend.position = "right",
plot.title = element_text(hjust=0.5)) +
scale_colour_manual(values = mycolors) +
labs(y = paste("Fitted", scalar), x = "Age", title=title2) + ylim(ylim1, ylim2)
}
return(ggarrange(plot1, plot2, ncol=2, common.legend=T, legend="right"))
}
plot_nodes_zeroed_hemi <- function(tract, df, scalar, ylim1, ylim2) {
df$nodeID <- as.numeric(df$nodeID)
if(str_detect(tract, "Forceps")) {
plot <- df %>% filter(tractID == tract) %>%
group_by(tract_node) %>%
ggplot(aes(x = age, y = est, group = tract_node)) +
geom_line(aes(colour = nodeID), size = 0.8, alpha = 0.6) +
theme_classic() +
theme(text=element_text(size=12),
legend.position = "none",
plot.title = element_text(hjust=0.5)) +
paletteer::scale_colour_paletteer_c("grDevices::RdYlBu", direction=-1) +
labs(x="", y="", title=tract) + ylim(ylim1, ylim2)
return(plot)
} else {
plot1 <- df %>% filter(tractID == tract) %>%
filter(hemi == "Left") %>%
group_by(tract_node) %>%
ggplot(aes(x = age, y = est, group = tract_node)) +
geom_line(aes(colour = nodeID), size = 0.8, alpha = 0.6) +
theme_classic() +
theme(text=element_text(size=12),
legend.position = "none",
plot.title = element_text(hjust=0.5)) +
paletteer::scale_colour_paletteer_c("grDevices::RdYlBu", direction=-1) +
labs(x="", y="", title=paste("Left", tract)) + ylim(ylim1, ylim2)
plot2 <- df %>% filter(tractID == tract) %>%
filter(hemi == "Right") %>%
group_by(tract_node) %>%
ggplot(aes(x = age, y = est, group = tract_node)) +
geom_line(aes(colour = nodeID), size = 0.8, alpha = 0.6) +
theme_classic() +
theme(text=element_text(size=12),
legend.position = "none",
plot.title = element_text(hjust=0.5)) +
paletteer::scale_colour_paletteer_c("grDevices::RdYlBu", direction=-1) +
labs(x="", y="", title=paste("Right", tract)) + ylim(ylim1, ylim2)
return(plot_grid(plot1, plot2, nrow=2))
}
}
plot_nodes_zeroed <- function(tract, df, scalar, ylim1, ylim2) {
df$nodeID <- as.numeric(df$nodeID)
plot <- df %>% filter(tractID == tract) %>%
group_by(tract_node) %>%
ggplot(aes(x = age, y = est, group = tract_node)) +
geom_line(aes(colour = nodeID), size = 0.8, alpha = 0.6) +
theme_classic() +
theme(text=element_text(size=12),
legend.position = "none",
plot.title = element_text(hjust=0.5)) +
paletteer::scale_colour_paletteer_c("grDevices::RdYlBu", direction=-1) +
labs(x="", y="", title=tract) + ylim(ylim1, ylim2)
return(plot)
}devtraj_zeroed_fa <- lapply(unique(nodewise_smoothfits_fa$tractID), plot_nodes_zeroed, df = nodewise_smoothfits_fa, scalar = "FA", -0.05, 0.03)
names(devtraj_zeroed_fa) <- unique(nodewise_smoothfits_fa$tractID)
devtraj_zeroed_fa_plots <- arrange_by_orientation(devtraj_zeroed_fa, "Fractional Anisotropy")
x.grob <- textGrob("Age (years)",
gp=gpar(col="black", fontsize=12))
y.grob <- textGrob("Zero-centered FA",
gp=gpar(col="black", fontsize=12), rot=90)
space.grob <- textGrob(" ",
gp=gpar(col="black", fontsize=12), rot=90)
grid.arrange(arrangeGrob(devtraj_zeroed_fa_plots, left = y.grob, bottom = x.grob, right = space.grob))
#ggsave("/Users/audluo/PennLINC/luowm_local/output/tract_profiles_testing/HCPD/devtraj_zeroed_fa.png", grid.arrange(arrangeGrob(devtraj_zeroed_fa_plots, left = y.grob, bottom = x.grob, right = space.grob)), width = 22, height = 24, units = "in")We see dissociable trajectories along each tract! - For A-P tracts, lower node IDs (blue; anterior segments) appear to have more steeply decreasing trajectories. Higher node IDs (red; posterior segments) appear to flatten out around mid-adolescence. Deeper tract regions (yellow) also have distinct trajectories that appear flatter than either end of the tract. - For R-L tracts, we also see distinct trajectories in deep white matter (yellow) compared to peripheral (red/blue), though forceps major seems to have flatter trajectories in left portions (blue) of its tract - For S-I tracts, deeper white matter (yellow and red for CST; yellow for the other tracts) have flatter developmental trajectories compared to that at tract ends (blue for CST; red and blue for other tracts)
devtraj_zeroed_md <- lapply(unique(nodewise_smoothfits_md$tractID), plot_nodes_zeroed, df = nodewise_smoothfits_md, scalar = "MD", -0.00005, 0.00006)
names(devtraj_zeroed_md) <- unique(nodewise_smoothfits_md$tractID)
devtraj_zeroed_md_plots <- arrange_by_orientation(devtraj_zeroed_md, "Mean Diffusivity")
x.grob <- textGrob("Age (years)",
gp=gpar(col="black", fontsize=12))
y.grob <- textGrob("Zero-centered MD",
gp=gpar(col="black", fontsize=12), rot=90)
space.grob <- textGrob(" ",
gp=gpar(col="black", fontsize=12), rot=90)
grid.arrange(arrangeGrob(devtraj_zeroed_md_plots, left = y.grob, bottom = x.grob, right = space.grob))
#ggsave("/Users/audluo/PennLINC/luowm_local/output/tract_profiles_testing/HCPD/devtraj_zeroed_md.png", grid.arrange(arrangeGrob(devtraj_zeroed_md_plots, left = y.grob, bottom = x.grob, right = space.grob)), width = 22, height = 20, units = "in")I also tried looking at developmental trajectories of nodes using
tractwise GAMs. So I fit my GAM for each tract:
gamm(MD ~ s(age, k = 3, fx = T, bs = 'cr') + s(nodeID, k = 24, fx = F, bs = 'cr') + ti(age, nodeID, k=c(3, 24), fx = F) + sex + mean_fd, random = list(subjectID=~1), data = df),
where I used bs=‘cr’ for the smooth spline basis to match the ti()
default spline basis. Then I used fitted_values to get the
fitted MD for each node at each age.
It’s definitely easier to see the differences in developmental trajectories with the zero-centered smooths above, but this is also kind of interesting… There’s a super stark difference between central (nodes 40-59) and peripheral (nodes 0-19 and 80-99) tract regions, where deeper white matter have super distinct MD profiles between tracts.
# Define the number of colors you want
nb.cols <- 7
mycolors <- colorRampPalette(paletteer::paletteer_d("LaCroixColoR::Pamplemousse"))(nb.cols)
AP_tract_ends_md <- plot_tract_ends("AP", fitted_md, "MD", mycolors, 0.00045, 0.0008)
nb.cols <- 5
mycolors <- colorRampPalette(paletteer::paletteer_d("LaCroixColoR::Pamplemousse"))(nb.cols)
AP_frontotemp_tract_ends_md <- plot_tract_ends("AP_frontal_temporal", fitted_md, "MD", mycolors, 0.00045, 0.0008)
SI_tract_ends_md <- plot_tract_ends("SI", fitted_md, "MD", mycolors, 0.00045, 0.0008)
nb.cols <- 10
mycolors <- colorRampPalette(paletteer::paletteer_d("LaCroixColoR::Pamplemousse"))(nb.cols)
CP_tract_ends_md <- plot_tract_ends("CP", fitted_md, "MD", mycolors, 0.00045, 0.0008)
ggarrange(AP_tract_ends_md, AP_frontotemp_tract_ends_md, SI_tract_ends_md, CP_tract_ends_md, nrow=2, ncol =2)
#ggsave("/Users/audluo/PennLINC/luowm_local/output/tract_profiles_testing/HCPD/fitted_val_tract_ends_md.png", ggarrange(AP_tract_ends_md, AP_frontotemp_tract_ends_md, SI_tract_ends_md, CP_tract_ends_md, nrow=2, ncol =2) + bgcolor("white"), width = 20, height = 10, units = "in")What do central vs. peripheral nodes look like for specific tracts?
scalar <- "MD"
df <- fitted_md
title1 <- "Peripheral"
title2 <- "Central"
ylim1 <- 0.00045
ylim2 <- 0.0008
plot_central_peripheral <- function(tract) {
plot1 <- df %>% filter(tractID == tract) %>%
filter(nodeID < 10 | nodeID > 90) %>% group_by(tract_node) %>%
ggplot(aes(x = age, y = fitted, group = tract_node)) +
geom_line(aes(colour = tractID), size = 0.9, alpha = 0.7) +
theme_classic() +
theme(text=element_text(size=16),
legend.position = "bottom",
plot.title = element_text(hjust=0.5)) +
scale_colour_manual(values = "#4CC3FFFF") +
labs(y = paste("Fitted", scalar), x = "Age", title=title1) + ylim(ylim1, ylim2)
plot2 <- df %>% filter(tractID == tract) %>%
filter(nodeID > 39 & nodeID < 60) %>% group_by(tract_node) %>%
ggplot(aes(x = age, y = fitted, group = tract_node)) +
geom_line(aes(colour = tractID), size = 0.9, alpha = 0.7) +
theme_classic() +
theme(text=element_text(size=16),
legend.position = "bottom",
plot.title = element_text(hjust=0.5)) +
scale_colour_manual(values = "#FF9932FF") +
labs(y = paste("Fitted", scalar), x = "Age", title=title2) + ylim(ylim1, ylim2)
return(plot_grid(plot1, plot2))
}
central_peripheral_plots <- lapply(c("Arcuate Fasciculus", "Cingulum Cingulate", "Forceps Major", "Forceps Minor", "Inferior Fronto-occipital Fasciculus", "Inferior Longitudinal Fasciculus", "Superior Longitudinal Fasciculus", "Uncinate Fasciculus", "Posterior Arcuate","Vertical Occipital Fasciculus"), plot_central_peripheral)
names(central_peripheral_plots) <- c("Arcuate Fasciculus", "Cingulum Cingulate", "Forceps Major", "Forceps Minor", "Inferior Fronto-occipital Fasciculus", "Inferior Longitudinal Fasciculus", "Superior Longitudinal Fasciculus", "Uncinate Fasciculus", "Posterior Arcuate","Vertical Occipital Fasciculus")
central_peripheral_specific <- ggarrange(plot_grid(plotlist=central_peripheral_plots[c(1:4)], ncol=2), plot_grid(plotlist=central_peripheral_plots[c(5:8)], ncol=2), plot_grid(plotlist=central_peripheral_plots[c(9:10)], ncol=2), heights = c(1,1,0.5), nrow=3)
central_peripheral_specific
#ggsave("/Users/audluo/PennLINC/luowm_local/output/tract_profiles_testing/HCPD/fitted_val_tract_ends_cp_md.png", central_peripheral_specific, width = 13, height = 20, units = "in")For some of the tracts, there’s kind of two chunks of developmental trajectories - could they correspond to each end of the tract?
scalar <- "MD"
df <- fitted_md
title1 <- "Anterior"
title2 <- "Posterior"
ylim1 <- 0.00045
ylim2 <- 0.0008
plot_anterior_posterior <- function(tract) {
plot1 <- df %>% filter(tractID == tract) %>%
filter(nodeID < 10) %>% group_by(tract_node) %>%
ggplot(aes(x = age, y = fitted, group = tract_node)) +
geom_line(aes(colour = tractID), size = 0.9, alpha = 0.7) +
theme_classic() +
theme(text=element_text(size=16),
legend.position = "bottom",
plot.title = element_text(hjust=0.5)) +
scale_colour_manual(values = "#4CC3FFFF") +
labs(y = paste("Fitted", scalar), x = "Age", title=title1) + ylim(ylim1, ylim2)
plot2 <- df %>% filter(tractID == tract) %>%
filter(nodeID > 89) %>% group_by(tract_node) %>%
ggplot(aes(x = age, y = fitted, group = tract_node)) +
geom_line(aes(colour = tractID), size = 0.9, alpha = 0.7) +
theme_classic() +
theme(text=element_text(size=16),
legend.position = "bottom",
plot.title = element_text(hjust=0.5)) +
scale_colour_manual(values = "#FF9932FF") +
labs(y = paste("Fitted", scalar), x = "Age", title=title2) + ylim(ylim1, ylim2)
return(plot_grid(plot1, plot2))
}
anterior_posterior_specific <- ggarrange(plot_anterior_posterior("Inferior Longitudinal Fasciculus"),
plot_anterior_posterior("Uncinate Fasciculus"),
plot_anterior_posterior("Inferior Fronto-occipital Fasciculus"), nrow=3)
anterior_posterior_specific
#ggsave("/Users/audluo/PennLINC/luowm_local/output/tract_profiles_testing/HCPD/fitted_val_tract_ends_ap_md.png", anterior_posterior_specific, width = 10, height = 15, units = "in")Other things to look at: - double check posterior-anterior and inferior-superior gradients - age of first developmental slowing - compute Euclidean distance to cortex; look at development vs. distance to cortex (instead of relative node position) - PCA on nodal age effect data?